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Summary 

This Quarterly Report deals with the three sets of topics summarized 
hereunder. For each set of topics the name of the PI responsible for that 
particular set is indicated. 

(l) E.C. Lorenzini, PI. 

Investigation of the propagation of longitudinal and transverse waves along 
the upper tether. Specifically the upper tether is modelled as three massive 

platforms connected by two perfectly elastic continua (tether segments). The 
tether attachment point to the station is assumed to vibrate both longitudinally 
and transversally at a given frequency. Longitudinal and transverse waves 

propagate along the tethers affecting the acceleration levels at the elevator and at 
the upper-platform. The displacement and acceleration frequency-response-func- 
tions at the elevator and at the upper-platform are computed for both longitudinal 
and transverse waves. 

An analysis to optimize the damping time of the longitudinal dampers is 

also carried out in order to select the “optimal” damper parameters. The 

analytical evaluation of the performance of tuned longitudinal dampers vs. detuned 


longitudinal dampers is also part of this analysis. 
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(2) G.E. Gullahorn, PI. 

Due to other demands on the time of one PI, reduced effort was spent on 
the tasks (A) for support of the Tether Applications Simulation Working Group and 
(B) for study of the use of the Shuttle primary Reaction Control System (RCS) 
thrusters for blowing away a recoiling broken tether. In support of the latter, 
a microcomputer system was acquired and set up. SLACK code has been 
transferred to this system from the VAX, and translated from VAX extended 
Fortran to the Fortran-77 standard used on the microcomputer. 

(3) R.D. Estes, PI. 

Most of the effort in the tether plasma physics study was devoted to 
software development in this period. A particle simulation code has been 
integrated into our Macintosh II computer system and will be utilized for studying 


the physics of hollow cathodes. 
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Figure Captions 

Figure 1. Schematic of the upper-tether. 

Figure 2(a)-2(f). Displacement and acceleration frequency-response- functions 

(FRF’s) at the elevator for longitudinal and transverse waves 
propagating in the upper-tether. 

Figure 3(a)-3(h). Same as Figures 2 except that the FRF’s are those of the 

upper-platform instead of the elevator. The longitudinal 
and transverse displacement attenuation-functions for waves 
crossing the elevator and reaching the upper-platform are 
shown in Figures 3(g) and 3(h). 

Figure 4. Schematic of a longitudinal damper and associated tether 

segment. 

Figures 5(a)-5(c). Longitudinal dynamic response vs. damper parameters and 

tether characteristics of tether segment 1 (between the 
lower-platforrn and the station). These figures are also 
representative of the dynamic response of tether segment 3 
(between the elevator and the upper-platform). 

Figures 6(a)-6(c). Same as in Figures 5 except that the dynamic response 

is that of tether segment 2 (between the station and the 
elevator). 

Figure 7. First “card” of hypercard front-end to simulation program. 

User defines the basic run parameters. 

Figure 8. Example of a card to define the input parameters for a given 

species of simulation particle. 
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1.0 INTRODUCTION 

This is Quarterly Report #13 submitted by the Smithsonian Astrophysical 
Observatory (SAO) under NASA/MSFC contract NAS8-36606, “Analytical Investi- 
gation of the Dynamics of Tethered Constellations in Earth Orbit (Phase II).” 
The PFs for this report are: Dr. Enrico C. Lorenzini for the analysis described in 
Sections 2.0-4.0, Dr. Gordon E. Gullahorn for Sections 5.0 — 7.0, and Dr. Robert D. 
Estes for Sections 8.0—10.0. This report covers the period from 1 April 1988 


through 30 June 1988. 
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2.0 TECHNICAL ACTIVITY DURING REPORTING PERIOD AND PRO- 
GRAM STATUS, E.C. LORENZINI, PI 

2.1 Wave Propagation In The Upper Tether Of TECS 

2.1.1 Introductory Remarks 

An important topic of investigation for the Tether Elevator /Crawler System 
(TECS) is the propagation of disturbances along the tethers. Specifically, the 
propagation along the upper tether is of particular relevance to the microgravity 
experiments on board the elevator [EL). The Space Station [SS), in fact, is a 
source of non-negligible disturbances at a wide range of frequencies. In general, 
the low frequencies (around 10"^ Hz) disturbances are associated with aerodynamic 
and orbital perturbations, the medium frequencies (lO'- Hz— 10 Hz) disturbances 
with the structural vibrations of the station, and the higher frequencies (>10 Hz) 
with rotating machinery and human activity on board the station. 

Since the tether is an elastic continuum the above mentioned disturbances 
travels from the SS to the upper-tether-end and back affecting the microgravity 
level on board the elevator. For a perfectly elastic tether with non-dissipative 
terminations these waves propagate indefinitely back and forth along the tether. 
In an actual case the tether has a small amount of material damping and is, 
moreover, embedded in a dissipative medium. The same applies to the 
terminations (platforms). 
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The amount of tether material damping has been estimated in a 
preliminary way. The tether is a complex non-isotropic continuum and the 
damping varies with the tether length according to a function which depends on 
the damping model adopted. Furthermore, material damping is significantly 
affected by temperature because the material properties of kevlar are sensitive to 
temperature variations. Current estimate of tether material damping ranges from 
1% to 5% of the critical damping for a tether length of 20 km [1,2]. The tether, 
therefore, has to be viewed as a low-damping elastic continuum. 

The small tether material damping affects primarily the propagation of 
longitudinal waves, while has an almost negligible effect on the transverse waves 
[l|. Transverse waves, on the other hand, are affected by the interaction with the 
surrounding atmosphere. The extent to which this interaction provides non- 
negligible damping of transverse waves along the tether has not been evaluated so 
far but it is quite reasonable that it is small given the low value of the tether 
transverse velocity. 

The first issue to be addressed, therefore, is the response of the upper 
tether, considered as a perfectly elastic continuum, when its attachment point 
to the station is oscillating in either the longitudinal or the transverse direction 
with a given frequency. 

If we consider the realistic case of small perturbations we can treat the 
longitudinal and the transverse waves independently as shown in the next 
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subsections. 


2.1.2 Transverse Waves 

The upper tether subsystem is schematically shown in Figure 1. The tether 
segment 1 of length £i connects the attachment point of the station to the elevator. 
The tether segment 2 of length £2 connects the elevator to the upper platform. 
The longitudinal coordinates along the undeformed tether segments 1 and 2 are 
called and z<i respectively; therefore 0 < z^ < ly and 0 < Z 2 < £ 2 - The 
attachment point at the station is perturbed in either the transverse or the 
longitudinal direction as indicated in Figure 1. In particular the transverse 
perturbation generates transverse waves propagating in the two tether segments. 
We call V'l and ^2 the lateral deformations of tether segments 1 and 2 produced by 
the transverse waves. 

After assuming that the platforms are point masses the transverse wave 
equations in tether segments 1 and 2, and the boundary conditions are given by 


dt^ 






dt^- 


H2 


dzl 


0 < < £1 

0 < Z2 < ^2 


( 1 . 1 ) 


(1.2) 


<h{0,t)= pre^‘ 


(1.3) 
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dz\ 


(^ 1,0 


r, - 




( 0,0 


r, = 


dz-i 

rjji{iut)= V»2(0,0 


dt^ 


( 0,0 


mi 


^02 


(^2,0 


T, = 


^^02 


(^2,0 


m2 


(1.4) 

(1.5) 

( 1 . 6 ) 


where we have expressed the external perturbation in equation (1.3) as a complex- 
variable function by utilizing Euler’s formula. The boundary condition of 
equations (1.4) states that the inertial force at the elevator (mi) is balanced by the 
resultant of the transverse forces generated by tether segments 1 and 2. Equation 
(1.6) is similar to equation (1.4) except that is referred to the upper platform (m 2 ) 
which is connected to tether segment 2 only. Equation (1.5) is an interface 
condition between the transverse displacement of tether segment 1 and 2 at the 
elevator. The solutions of equations (l.l) and (1.2) have the form 


01 = i?i(2i)e-‘ 

(2.1) 

02 = -^ 2 ( 22 ) 6 '“* 

(2.2) 

After substituting equations (2) into equations (l.l) and (1.2) we 
differential equations in /Zi( 2 i) and whose solutions are 

obtain standard 

Rj{zj) = ] j — 1,2 

( 3 ) 
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where Dy’s are complex numbers which can be expressed as follows 

Dy = dj + idj J = 1,2 (4) 

and 

kj = uj/c^j i = 1,2 (5) 

We define 

<f>j = J = 1,2 (6) 

and after substituting equations (3) into equations (2) and subsequently into the 
boundary conditions (1.3) — (1.6) we obtain 


di = pt 

^djcos^i + dism 4>i^Tiki + (^d2COS 4>2 d2s\n 4>2^Tik2 = - u^mid2 

di cos <j)i — d[ sin <j>]^ — d^ 

^^2 cos ^>2 + sin (/» 2 ^ T’ 2^2 = — w^m 2 ^ ^2 (^2 ~ ^2 sin 02 ) 


( 7 ) 


By solving equations (7) for d^, dj, d 2 , and d 2 we can explicit the shape functions 
Rj{zj) given by equations (3). In particular we are interested in the values of the 
function /?i(2i) and R{z 2 ) at the coordinates Z[ = ti and Z 2 = ^2 respectively. 
Ri{ii) and i?2(^2), in provide the transverse displacements at the elevator and 
at the upper platform respectively for a transverse displacement perturbation of 
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the tether attachment point to the station. 

After several algebraic manipulations we obtain the displacements at the 
elevator and at the upper platform as follows: 


where 


1 

and 


— Pt^i 

(7.1) 

^2(^2) Pt ’^2 = 

(7.2) 

, . , , cuoo sin d>i 

cos (pi - cjUi sm <pi + T — — — 

-I 

(8.1) 

cos q>2 — u}02 sin (p2 


\ 6 t\ — 1 cos (f>2 - u}Q2 sin (j>2 \ 

(8.2) 

IX2I = |Xx||^r| 

(8.3) 

ay = rrijC^j/Tj j = 1,2 

(9.1) 

4 >j = w£y/c^y j = 1,2 

(9.2) 

T = T 2 C^J{TiC^^) 

(9.3) 


I Xx I and | X 2 | are the transverse-displacement frequency-response-functions [FRF) 
at the elevator and at the upper platform. From equations (7) we infer that | | 

provides the attenuation of a transverse wave passing through the elevator and 
tether segment 2. It is comforting to notice that for T 2 = 0 (no wave propagation 
in tether 2) equation (8.1) reduces to the form shown in reference [2] which has 
been derived for a single tether seg^uent. Since the wave equations are linear we 
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can easily obtain the acceleration FRF^s Xi and X 2 as follows [l] 


X, 


4>j I Xy 


j = 1,2 


( 10 ) 


2.1.3 Longitudinal Waves 

By inspecting equations (l) we see that the transverse wave equations and 
their boundary conditions can be transformed into longitudinal wave equations and 
associated boundary conditions by making the following substitutions: (a) the 

transverse displacement iij is substituted by the longitudinal displacement qy; (b) 
the transverse wave velocity by the longitudinal wave velocity c,,y; (c) the 

tensions Tj by (EA)j; and (d) the transverse displacement amplitude pf by the 
longitudinal displacement amplitude pi. After the transformation the longitudinal 
wave equations are obtained as follows: 


dt^ 



0 < Zj < l.j, j = 1,2 


( 11 . 1 ) 


qi(0,<) = 


( 11 . 2 ) 


drn 

dzi 


i^ut) 


{EA), 


dri2 

dZ2 


(0,t) 




d^V2 

dt^ 


( 0,0 


mi 


(11.3) 


»7i(0,0 = 


(11.4) 
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dm 

dz2 


(^2,0 


{EA)2 


d^m 


(^25 0 


m2 


(11.5) 


where the longitudinal wave velocities are given by: 


S = \J(EA)j/tij i = 1,2 


( 12 ) 


The displacement FRF's | | at the elevator and | A2 | at the upper-platform, and 

the attenuation function | di j for the longitudinal waves are also obtained from 
equations (8) by performing the above mentioned transformation. We have 
therefore 



where 


cos ^71 - ubi sin 7i -h ^ 


u>b 2 sin 7i 


cos 72 " ^^2 sin 72 
\h\ - I cos 72 - 0162 sin 72 




A 2 — Aj 6j 


bj = mjC„j/{EA)j j = 1,2 

Vj = ^^j/% j = 1.2 

c ^ [EA)2C,,^ 

[EA\Cr,2 


(13.1) 

(13.2) 

(13.3) 

(14.1) 

(14.2) 

(14.3) 


The longitudinal displacements at the elevator and S2 at the upper platform are 
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therefore given by 


SM = PlAi (15.1) 

52 (^ 2 ) = PiAi^t (15.2) 


The acceleration FRF’s | A‘i | and | A 2 | are immediately obtained from equations 
( 13 . 1 ) and (13. 3) as 




3 = 1,2 


(16) 


2.1.4 Numerical Results 

The longitudinal and transverse- wave frequency response functions at the 
elevator and at the upper platform have been evaluated for the typical design 
parameters of the upper tether of TECS. Specifically (see Figure l) we have: 
{EA)i = [EA )2 = 61645 N, Ti = 394.5 N, = 375.5 N, and ^ 

4.9x10“^ kg/m. Consequently the transverse-wave velocities in tether segment 1 
and 2 (of the upper tether) are 


= 283.74 m/s 
= 276.83 m/s 


(17.1) 

(17.2) 
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and the longitudinal-wave velocities 

c^i = c ^2 = 3546.92 m/s (18) 

The elevator displacement FRF^s | Xj | (transverse waves) and | Aj | (longitudinal 
waves) are shown in Figures 2(a) and 2(b). Similarly the elevator acceleration 
FRF's I Xi 1 and | Ai | for the transverse and longitudinal waves are shown in 
Figures 2(c) and 2(d) respectively. The last two figures have been zoomed-in 
for frequencies less than 10 Hz in Figures 2(e) and 2(f). In general the tether 
segment behaves as the superposition of a single DOF low-pass filter and a 
transmission line for both longitudinal and transverse disturbances. In this 
analysis the effect of damping, both internal and external, has been neglected. In 
an actual case the resonance peaks will be limited by the damping. The lower 
envelopes of the displacement FRF's correspond to the response of a single DOF 
low-pass filter while the additional spectral lines are related to the transverse and 
longitudinal waves propagating in the tether. 

Figures 3(a)-3(f) are like Figures 2(a)-2(f) except that the former set is 
referred to the upper platform. By comparing these two sets of figures it is clear 
that the elevator acts as a wave attenuator. The noise at the upper platform is a 
few orders of magnitude smaller than at the elevator. Specifically Figures 3(g) 
and 3(h) show the displacement attenuation functions | | and | | for 

transverse and longitudinal wa’’- respectively. The acceleration attenuation 
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functions 0^ \ and 9‘1 are readily obtained as follows: 


df 1 

= (<^2/0l)' 1 Ot I 

(19.1) 

€\ 

= {li/'li? 1 1 

(19.2) 


Figures 3(g) and 3(h), therefore, provide also the acceleration attenuation functions 
when their magnitudes are scaled by 





{li/ll? 



81 


( 20 . 1 ) 


( 20 . 2 ) 


The response at the upper platform stresses the previous conclusion: the upper 

tether acts as the superposition of a low-pass filter and a transmission line. 
The elevator, however, attenuates the transverse and the longitudinal waves 
propagating from tether segment 1 (below the elevator) to tether segment 2 (above 
the elevator). 


I 



FRF Magnitude (Log) Magnitude (Log) 
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Figures 3(a) and 3(b) 
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2.1.5 References To Section 2.1 

1. Xiaohua He, and J.D. Powell, “Tether Damping in Space,” to appear in 

t 

Proceedings of the Second International Conference on Tethers in Space, 
Venice, Italy, 4-8 October 1987. 

2. G.E. Gullahorn and R.G. Hohlfeld, “Tether as a Dynamic Transmission Line,” 
to appear in Proceedings of the Second International Conference on Tethers in 
Space, Venice, Italy, 4-8 October 1987. 
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2.2 Longitudinal Dampers 
2.2.1 Introductory Remarks 

The damping of longitudinal (along the tether segments) oscillations is 
essential to provide acceleration levels on board the elevator and the station 
suitable for microgravity experiments [l]. Specifically once the deployment- 
induced-oscillations have been dissipated the libration and lateral oscillation control 
systems can be switched off without impairing the microgravity experiments [l]. 
In fact the steady-state external perturbations are not strong enough to excite 
appreciably the libration and the lateral oscillations. On the contrary, thermal 
perturbations have a non-negligible effect upon the longitudinal oscillations. 
Everytime the system crosses the terminator, the tethers’ temperatures vary 
abruptly (almost impulsively) and longitudinal oscillations are excited twice per 
orbit. Longitudinal oscillations, therefore, must be continuously controlled if low 
acceleration levels are desired. 

In reference [2] we have described the three longitudinal dampers which 
are in series with the three tether segments connecting massive platforms. In 
reference [2] we also selected each damper parameters (stiffness and damping 
coefficient) based on considerations arising from physical intuition. Each damper 
has been tuned to the frequency of the natural bobbing frequency of the associated 
tether segment. We also adopted a damping coefficient which provides a damping 
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ratio equal to 0.9 for an ideal damper which is decoupled from the associated 
tether (i.e. the ideal damper is directly connected to the platform). This selection 
of parameters has indeed provided effective damping of tether bobbing oscillations. 
Nevertheless, since the performance of the longitudinal dampers is essential to 
guarantee a microgravity environment on board the tethered system, we want to 
verify more rigorously the validity of our selection of the damper parameters. We 
also want to assess if a detuned damper (i.e. damper stiffness different from the 
tether stiffness) can possibly perform better than a tuned longitudinal damper. 


2.2.2 Optimization Of Longitudinal Damper Parameters 

A relatively accurate mathematical model, amenable to an analytical 
solution, for each damper and the associated vibrating tether segment is as follows: 
each damper is assumed massless and only two degrees of freedom (DOF) are 
taken into account, i.e. the damper stretch and the elastic stretch of the associated 
tether segment (see Figure 4). The assumption of massless damper is actually 
quite accurate for an active control system. As a result of high feedback gains, 
the active damper’s response exhibits a negligible phase delay and therefore the 
effects of the damper’s (spool’s) inertia are masked. The two-DOF-assumption 
provides a realistic picture of the actual behavior of the system for small values of 
lateral oscillations since the coupling between the lateral and the longitudinal 
oscillations decreases quadratically w’‘h the lateral oscillation amplitude. This case 
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of “quasi "-aligned tethered system is actually the most interesting because is 
representative of the system steady-state configuration. 

With reference to Figure 4 the force-balance equations for the longitudinal 
vibrations of one tether segment with damper are as follows: 


— F = Fg (21.1) 

F = F,i (21.2) 


where the equivalent mass rriq is equal to mj • m‘i/[rni -f m 2 ), L is the geometrical 
distance between two adjacent platforms with masses mj and mo, Fg is an 
arbitrary external force. The forces F (tether) and F,i (damper) are given by: 


F = - = - kit (22.1) 

= -[k,il,i + bla) ( 22 . 2 ) 


where l„ is the unstretched tether length, l,i is the damper stretch, k is the tether 
stiffness, k,i and b are the stiffness and damping coefficient of the damper, and It is 
the elastic stretch of the tether. 

After substitution of equations (22) into equations (21) and after some 
algebraic manipulations, we obtain the equations of motion as follows: 
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£ + cj|£ - cj^d = Fe/ttiq (23.1) 

Id + 2X,/£,/ - X£ = 0 (23.2) 


where t = L - to the total stretch of the vibrating system. We have also 
defined 


= k/rriQ 

Xd = kd/h (24) 

X = k/ b 

Equations (23) describe the motion of a 3'‘*-order vibrating system which is 
actually amenable to an analytic solution. 

By Laplace transforming equations (23), with null initial conditions, we 

obtain 


s 2[£| + wj(£| - w, “[£,,] = \Fe\/ mq 


(25.1) 


sl£,,]+2X,/l£,/]-X[£]= 0 


(25.2) 


where [ ] denotes the Laplace transformation. After solving equation (25.2) for 
[£,l] and substituting into equation (25.1) we obtain 
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m 


X 

s + 2X,i 


[^1 


(26.1) 


[£]{s2 + h? - w,^X/(s + 2X,/)} = \FE]/mQ 


(26.2) 


which after substitution into equation (25.1) provides the characteristic equation of 
the system 


+ 2X,,s^ + ui’(2X,, -X) = 0 


(27) 


By using the transformation s = z - 2X,//3 equation (27) reduces to the canonical 
form 


+ pz + q - 0 


(28) 


where 


P = 

q = u;2|^|x,/-x^- A X,3 (29) 


The discriminant of the cubic equation (28) is 



+ 1 ? 


A 


3 


(30) 
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For an oscillatory response of the system we must have > 0 which implies that 
a pair of roots are complex conjugate and one root is real. For Zl < 0 the 
roots are all real. Real roots of equation (28) implies real roots of equation (29). 
By applying the Routh criterion to equation (29), however, we discover that at 
least one of the three real roots is positive. Consequently the system is unstable 
for A 0 and we limit our parametric analysis to the case Zi ^ 0. 

Instead of using Cardano’s formulae (which are quite cumbersome) for 
computing the roots of the cubic equation (27), we compute the real root by 
means of a numerical root finding routine and then we deflate equation (27). 
After deflation, equation (27) becomes 

(s — Si) [s^ + (S 2 + 2X,i)s + (Sj + 2X,/)5j + c]= 0 (3l) 

Since equation (31) has real coefficients, the roots S2 S3 of the 2“‘*-order 

equation between square brackets are complex conjugate. After defining 

6*2 — “1“ 

S3 = /? - iw (32) 

we can express the characteristic equation (27) as 


D{s) = (s - S|) [ . - (/9 + iw)] [s - (/? - ia>)] 


(33) 
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From the solution of the quadratic equation between square brackets in equation 
(31) we have 


(3 = -{si + 2X,i]/2 

r 1^2 

u = -|(si + 2X,/)^ - 4[{si + 2X,/)si + (34) 

We can express [£] and [£,i] in terms of transfer functions by using equations 

(26) 

[£] = (35.1) 

l^,/l = G,{s)\FE]/mc^ (35.2) 

The transfer functions are given by 

6'i(5) = {s^2Xu)/D{s) (36.1) 

G,(s) = X/D{s) (36.2) 

where D{s) is the characteristic equation. 


The expressions of £ and £,/ as a function of time [i.e. the solution of the 
3“horder differential equations (23)] are obtained by Laplace reverse-transforming 
equations (35). A case of part. -ular theoretical and practical interest is the 
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response of the system to an impulse. If we assume that Fg is an impulsive 
function of strength I we have that 


[Fjj] = / = constant 


(37) 


and equations (35) can be Laplace reverse-transformed by a lengthy but well 
known technique. Gi( 5 ) and G(s) are first reduced to a sum of simple fractions as 
follows 


Gi{s) 

<^2(5) 


5 - Si £^2(5) D‘i[s) 

= 4. ^2 

s - Si D^is) D<i{s) 


(38.1) 

(38.2) 


where £> 2 (^) = ~ (/^ + *^)] [^ “ (/? - »w)] and the coefficients Ai, Bi, Ai, 

B[, and B 2 are functions of X,/,X and u„. 

The variable of greatest interest to us is actually the elastic stretch 

£1 = £-£,1 (39) 


Hence, by taking into account that the Laplace transform is a linear operator, we 
can define a third transfer function (^ 3 ( 5 ) as follows: 
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[^t] = Gz{s)I/mq (40.l) 

where 

<^ 3 ( 5 ) = G,(5)-G2(s) (40.2) 

Subtracting equation (38.1) from equation (38.2) we obtain the fraction formulation 
of ^ 3 ( 5 ) and we can finally Laplace reverse-transform equation (40.1). We give 
only the final result, without providing details on the reverse-transformation for 
the sake of brevity. 

The elastic stretch as a function of time for null initial conditions and 
impulsive external force is given by 

[d cos(wt) + B sin(cui)]| (4l) 

where /? and cu have already been defined and 


2X,/ — X Sj 
2 •> 

3sj -|- + 4X,/Sj 


B = 



(42) 


Equation (41) gives the response >f the elastic stretch as a function of the 
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parameters of the 3‘^‘^-order vibrating system X, X^, and 

Since the acceleration measured on any platform of the tethered system is 
directly proportional to the tether elastic stretch, the time history of the latter 
quantity is of particular significance to the microgravity levels of our tethered 
system. In an actual case a tethered system is excited twice per orbit by a 
“quasi” impulsive thermal perturbation which deposes energy primarily into the 
longitudinal oscillations. Our goal is therefore to define those values of parameters 
X, X,i and u}„ which provides the smallest and shortest fluctuation of the 
acceleration. 

Several parameters can be used as indicators of the effectiveness of the 
damper. We adopt the settle time f.v and the maximum elastic stretch per unit 
initial velocity Et [in equation (40.1) I/mQ is equal to the initial velocity] as 
primary indicators of the damper’s effectiveness. We define the settle time as the 
time taken by the system to reduce the oscillation amplitude to 10% of its 
maximum value after the application of the external impulse. We also define the 
rise time as the time to reach the maximum amplitude of the response and the 
relaxation time tp as the time the system takes to reduce the amplitude from its 
maximum value to 33% of its maximum amplitude (i.e. ^(/e). 

The rise time and the maximum elastic stretch have been computed for 
several values of X and X,/ by bracketing equation (40.1) around its absolute 
maximum and then using a root-fine.. ^ig routine to locate the maximum value with 
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high accuracy. The computation of the settle time and relaxation time require 
knowledge of the envelope of tf Such an envelope is given by: 




mq 



-j- M 


) 


(43) 


( \'''‘ 

where M = ( A^+B^j . The relaxation time is then computed by solving 
numerically the following equation 

I |e.n(i/J+‘M) + I M = Ijt (44) 

where is the previously computed rise time. The settle time <v is also computed 
by equation (44) after replacing the factor e with 10. We want to point out that 
in general + tji. 

Figures 5(a)-5(c) depict the settle time, the maximum elastic stretch per 
unit initial velocity, and the pseudo frequency vs. X and X,/ for = 0.0242 

rad/sec (/, = 3.844 x 10“^ Hz). This value of u>„ is the natural bobbing frequency 
of the first tether segment (i.e. between the lower mass and the station). 

The settle time, in Figure 5(a), depends almost linearly upon the damper 
parameter X,/. The dashed line connects the points with X=X,/ which are 
representative of a damper tuned to the bobbing frequency of the associated tether 
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segment. The shortest settle time for a tuned damper is obtained for 
X = Xfi ~ 0.0135 sec'^ which is actually the value adopted in reference [6] for the 
longitudinal damper of tether 1. The settle time for a detuned damper can be 
shorter than that for a tuned damper as shown by all the points below the dashed 
line. From inspection of Figure 5(b) we notice, however, that some of the above 
mentioned points (for a detuned damper) correspond to maximum elastic stretches 
much larger than those associated with a tuned damper. There are only small 
areas in the X-X,i plane which correspond to a short settle time and a small 
elastic stretch for a detuned damper. One of such areas is for those values of X 
and X,i along the lowest points of Figures 5(a) and 5(b). Below those points the 
response is unstable and a certain margin should be provided to account for 
possible fluctuations of X and X,^ 

Figure 5(c) shows the pseudo-frequency a; as a function of X and X,i for 
(jjf, = 0.0242 rad/sec. It is worth noticing the reduction of the pseudo-frequency 
for certain values of X and X,/ of a detuned damper. 

Figures 6(a)-6(c) depict the settle time, the elastic stretch, and the pseudo 
frequency respectively vs. X and X,/ for = 0.117 rad/sec (/„ = 1.868x 10“^ Hz) 
which is the natural bobbing frequency of tether segment 2 (between the station 
and a 1-km-off elevator). The trends of these plots are similar to the previous set 
of plots except for the ranges of the parameters X and X,/ which are larger than 
before. There may be a scaling factor related to w,, even if it is not immediately 
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evident from the inspection of the equations. 

This time the shortest settle time for a tuned damper is even closer to the 
best results obtained by adopting a detuned damper. The detuned damper is also 
affected, for certain values of the parameters, by the problem of the amplification 
of the elastic stretch. We prove again that the values adopted in reference [6] for 
the longitudinal damper of tether 2, X = X,/ = 0.065 sec'^, provides the shortest 
settle time for a tuned damper. The figures for tether segment 3 are not shown 
because they are quite similar to Figures 5(a)-5(c) of tether segment 1. The 
natural frequency of tether segment 3 is - 0.0247 rad/sec (/„ = 3.934x 10“^ 
Hz) which is very close to the natural frequency of tether segment 1. In 
conclusion, because of the moderate advantages of a detuned damper over a tuned 
damper we have decided to adopt tuned longitudinal dampers, as in reference [l], 
with the following values of parameters 

Xi = X,,j = 0.0135 sec-^ 

Xo = X,/2 = 0.065 sec-^ (45) 

X 3 = X ,,3 = 0.0137 sec-^ 

The indexes in equations (45) refer to the three tether segments respectively. 
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2.2.3 References To Section 2.2 

1. E.C. Lorenzini, M. Cosmo, S. Vetrella and A. Moccia, “Acceleration Levels on 
Board the Space Station and a Tethered Elevator for Micro and Variable- 
Gravity Applications,” to appear in Proceedings of 2“^-International Conference 
on Tethers in Space,” Venice, Italy, 4-8 October 1987. 

2. E.C. Lorenzini et al., “Analytical Investigation of the Dynamics of Tethered 
Constellations in Earth Orbit (Phase II),” SAO Quarterly Report ^2 under 
contract NASA/MSFC NAS8-36606. 

2.3 Concluding Remarks 

A perturbation originated from the space station produces frequency 
responses at the elevator and at the upper-platform which contain frequency 
components related to the longitudinal and transverse waves propagating in the 
upper-tether. These frequency responses can be viewed as the superposition of the 
responses of a single-degree-of- freedom low-pass filter and of a transmission line. 
This result is not discouraging because the higher frequency longitudinal waves 
could be abated by the longitudinal dampers which are tuned to the low frequency 
bobbing frequencies. The abatement of transverse waves, on the other hand, 
requires a low-frequency transverse wave attenuator. By analyzing the results of 
our analysis we notice that the elevator attenuates strongly both the longitudinal 
and transverse waves travelling fror>" the tether segment below the elevator to the 
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one above. This result provides an indication of how to design a passive wave- 
attenuator for transverse waves which will be investigated during the next 
reporting period. 

The wave propagation analysis stresses once more the importance of the 
longitudinal dampers not only as devices for damping out the bobbing oscillations 
but also for attenuating the longitudinal waves. We have, therefore, carried out an 
analysis to optimize the damping performance of the longitudinal dampers. The 
results of this analysis shows that detuned dampers may be actually even more 
effective than tuned dampers in abating longitudinal oscillations. The advantages 
are, however, quite limited. Our analysis has confirmed that the damper 
parameters, that we had selected for the dampers, are “optimal” with respect to 
the damping time for tuned dampers. Since the bobbing frequencies are the lowest 
frequencies of the longitudinal wave frequency-response-functions the longitudinal 
dampers are also expected to act as low-pass filters. These issues related to the 
attenuation of longitudinal and transverse waves will be investigated in the next 


reporting period. 
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3.0 PROBLEMS ENCOUNTERED DURING REPORT PERIOD, E.C. LOREN- 
ZINI, PI 

None 

4.0 ACTIVITY PLANNED FOR NEXT REPORTING PERIOD, E.C. LOREN- 
ZINI, PI 

In the next reporting period we will analyze the techniques for attenuating 
the longitudinal and transverse wave propagating from the space station to the 
upper-platform along the upper tether. 

We will also work on the implementation of the rotational dynamics of the 
station and the elevator into MASTER20 computer code. 
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5.0 TECHNICAL ACTIVITY DURING REPORTING PERIOD AND PRO- 
GRAM STATUS, G. GULLAHORN PI. 

A variety of other demands (response to time critical demands of the TSSl 
mission as PI for an experiment on board, startup of a new grant, and internal 
administrative duties) reduced the effort the PI could devote to these tasks. 
Primary activity revolved around acquisition and setup of a new microcomputer 
system and transporting the SLACK code to that system. It is expected that 
greater effort can be devoted in the next reporting period. 


5.1 Tether Applications Simulation Working Group Support 
No significant activity in this reporting period. 
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5.2 Tether Aerodynamic Effect Of RCS Thruster Plume 

The activity during the reporting period involved acquisition and setup of a 
microcomputer system and transfer of Fortran code for the SLACK program to 
the new system. 

5.2.1 Computer Acquisition 

A Northgate microcomputer system with a 12 MHz 80286 central processor 
and 8 MHz 80287 numeric coprocessor was acquired. The system includes 1 MB 
of system memory, a 60 MB hard disk, and a monochrome graphics display. A 
variety of software was either transferred from a home computer system or may 
be used on both systems under manufacturer’s licensing agreements: compilers 

for Fortran, Pascal and C; wordprocessors; math libraries. Some further utility 
software was purchased: Pascal debugger; multitasking and task switching 

utilities; an 8087 library for the Fortran. 

It should be noted that an independently clocked 12 MHz 80287 (Microway) 
can be obtained at reasonable cost which would increase floating point throughput 
by up to 50%. We have also requested institutional (overhead) funds to upgrade 
the basic system to an 80386/80387 system, which should about double the system 
throughput initially, and with proper compilers (taking advantage of instructions 
and registers unique to the 386/387) increase floating point throughput by factors 
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of 4 to 10. 

5.2.2 Transport Of SLACK Code To Microcomputer 

Physical transfer of the source code for the SLACK program involved little 
difficulty; indeed, it was already available on diskette, having been transferred via 
modem for distribution to those requesting the source. 

Compilation and linking of the program proved not so simple. When 
written on the VAX full use had been made of a number of convenient extensions 
to the Fortran-77 standard provided in the VAX compiler (and one non-standard 
feature which though convenient is not simply an extension and which can lead to 
subtle errors). The compiler available (Microsoft) on the microsystem has few of 
these extensions, and a review of other microcomputer Fortrans showed that none 
of them supported the full set of extensions used. Hence a program FILTER was 
written (about 400 lines of Turbo Pascal) which accepts as input a VAX Fortran 
subroutine and outputs a version much closer to the F77 standard, flagging some 
constructs for manual modification. FILTER does the following 

• Removes in-line comments to a separate line. VAX Fortran allows a 
comment to be placed on the same line with code, set off after an 
exclamation point. F77 requires all comments to be on a separate line 
with a C in the first column. 

• Converts continuation lines starting with TAB to begin with 5 spaces. 

• Convert continuation line with continuation character ‘0’ to continuation 
character T’. Microsoft does not allow ‘O’. 
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• Convert initial TAB to 7 spaces. Microsoft would accept these lines, 
but this is done for consistency with the above. 

• Convert DO / ENDDO pairs to statement label format. VAX Fortran 
allows an extension in which the end of the range of a DO statement 
is delimited by an ‘ENDDO’, while F77 requires a labeled statement 
(referred to in the DO statement itself) as the delimiter. A statement 
label is provided (a five digit number), inserted in the DO statement, 
and attached to a created CONTINUE statement at the end of the range. 

• Flag DO WHILE statements for manual modification. This useful 
extension repeats a loop while some condition is met (e.g., DO WHILE ( X 
> 1.0 ) ). This construct is replaced by a conditional branch around the 
loop and an unconditional branch at the end. It proved easier to simply 
fiag these and modify them by hand. 

• Flag DATA and COMMONS statements for possible manual (a) transfer to 

BLOCK DATA (b) change of the commons name if same as subroutine 
name. VAX Fortran allows a labeled commons to have the same name as 
a subroutine; F77 requires all program units (including commons and 
subroutines) to have distinct names. VAX Fortran allows initializing of 
variables in commons with a DATA statement. F77 requires that this be 
done in a separate BLOCK DATA subprogram, and on reflection the VAX 
permissibility, though convenient, can lead to subtle bugs: if data 

statements are used to initialize the same variable in more than one 
subroutine, tests show that the linker does not pick up this inconsistency, 
and the value assigned depends on the order in which routines are linked. 
Good programming practice would: initialize variables local to a 

subroutine with data statements and not put them in commons, as this 
is not needed in the F77 standard (the habit of putting local variables 
in commons when their value must be retained between subroutine calls 
was developed using a non-standard, stack-oriented compiler); use 
parameter statement initialization where values need not be modified; use 
block data or other distinct initialization (as in a startup routine) where 
variables must be referred to by several routines. 


The SLACK code was run through the filter program and manual 
modification was made as needed; some cleanup could still be done on unneeded 
commons statements, as noted above. An additional problem arose in a subroutine 
REFINE which refines solutions to a single non-linear equation. This contains 
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convergence criteria in exceptional cases (which arise fairly often in SLACK) which 
depend on the details of real arithmetic. Two parameters (the maximum real 
number and a constant related to maximum precision) had to be modified since 
the VAX 8-byte real numbers are not the IEEE standard 8-byte reals used in 
virtually all current microcomputer compilers, whether for 8086 or 68000 series 
computers. The maximum real number for the IEEE format is readily available: 
about 1.797 x 10^°®. The precision factor 6, defined as the smallest number for 
which x 7 ^ (l+6)x can be assured (i.e. is true for all x), is not so easily found. It 
may be estimated from the length of the fraction part of the number, but its 
precise value depends on the details of rounding which are not given in typical 
compiler documentation. Estimates put 6 at 2'^^ or 2.22 x lO'^*^. A program was 
written to locate an “individual” S by logarithmic bisection for randomly chosen 
numbers x in an interval; the maximum individual S then forms the universal 6 
which will insure that the criterion is satisfied. Several runs for a variety of 

t 

ranges lead to 6 = 1.12 x 10'^^. 

A working version has been produced. It’s results need to be compared to 
those from the VAX and any discrepancies accounted for; some may be expected 
due to the difference in precision of the arithmetic (the IEEE standard has greater 
range than the VAX and about a digit less precision). The VAX version had been 
deleted due to disk space limitations during an extended period in which SLACK 
was not being used. It must either be restored from tape or recreated from source 
code, and a controlled set of comp :isons run. 



Page 49 


It must be noted that experience with the Microsoft compiler has been less 
than satisfactory. Fairly innocuous errors cause the computer to hang, requiring a 
coldboot. And although a working version has been created, the difference 
between that version and versions which fail to run (indeed hang the system) is 
unclear, seeming to reside in the options chosen on compiling (options chosen, 
indeed, to provide diagnostics on failure). Even when failing in a nondestructive 
way, the diagnostics provided are not very informative. Running with the 
CodeView debugger provided, though it is in many ways a very attractive tool, 
sometimes seems to generate failures or fail in different ways than without the 
debugger. Although all these problems arise occasionally with compilers on larger 
systems, the fragility of the Microsoft system supports the various word of mouth 
reports that it is not the preferred microcomputer compiler. (To be fair, most 
benchmarks show the Microsoft compiler producing moderately to substantially 
faster code than its rivals.) (Microsoft was initially chosen because of its 
attractive street price, about half that of other comparable compilers, a significant 
factor for personal acquisition.) If the problems experienced do not quickly 
resolve themselves, another compiler will be obtained; the two most generally 
favored are Lahey and Ryan McFarland, each with different advantages. 
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5.2.3 Modification Of SLACK Code 

Only slight effort was devoted to modifications of the SLACK code to reflect 
the more general drag in the trajectory calculation. Most effort was spent in 
getting the existing code working on the microcomputer; modifications will be 
made in that setting. Several other factors require analysis and inclusion in the 
code; periodic updating of the drag forces; specification of the thruster positions 
and firing strategy; calculation of the tether area presented, and tether 
configuration for lift; computing the aerodynamic forces due to the RCS plumes. 

6.0 PROBLEMS ENCOUNTERED DURING REPORTING PERIOD, G. GUL- 
LAHORN PI. 

None 
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7.0 ACTIVITY PLANNED FOR NEXT REPORTING PERIOD, G. GUL- 
LAHORN PI. 

7.1 Tether Applications Simulation Working Group Support 

Preparing a final report on this task will be the first priority during the 
initial period of the next quarter. A major effort for the report will be preparing 
a uniform set of comparison plots for the various cases and simulators; this effort 
is complicated by the fact that the Imagen laser printer which has been used for 
plotting does not have sufficient memory to process some of the plots needed; a 
second printer exists with (presumably) adequate memory, but different scalings 
and fonts will require modification of the plotting schemes which had been 
developed. An alternative is the use of a microcomputer based graphics system, 
which would typically allow more convenient placement of labels, etc. A number 
of such systems are being evaluated to see if they will support the number of lines 
and data points needed for our purposes; a high resolution dot matrix printer is 
available. Work on this task, other than reporting, is essentially complete. Some 
few further efforts may be made as suggested in Quarterly 11. 
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7.2 Tether Aerodynamic Effect Of RCS Thruster Plume 

The topic of aerodynamic forces on the tether due to the RCS plumes 
should be given further consideration, in particular to find an indication of the 
orthogonal (lift) forces. Experimental results may now exist. The relationship of 
the orthogonal force to the angle between the tether and the flow needs to be 
better understood. 

Modification of the SLACK code will be continued. It is probable that 
a preliminary version, at least, can be completed in the next quarter. 

The port of the SLACK code to the microcomputer will be studied and 
completed. The VAX version of SLACK will be resurrected, and comparison runs 
made. If the current compiler cannot be made to perform in a consistent and 
robust fashion, a new compiler may be acquired. 
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8.0 TECHNICAL ACTIVITY DURING REPORTING PERIOD AND PRO- 
GRAM STATUS, R.D. ESTES, PI 

Most of our effort in the study of the plasma physics of hollow cathodes for 
use with tethered satellite systems in this reporting period was devoted to software 
development. This was carried out on one of the Macintosh II computer systems 
recently purchased by the group, partially out of funds from this contract. The 
development system used was Lightspeed C, an integrated system consisting of an 
editor, C compiler, and linker that can generate “stand-alone” Macintosh 
applications which can be run as regular Macintosh programs on any Macintosh 
computer. The system has full support for all of the features of the Macintosh 
computer that make it “user friendly”, including menus, event-driven programs, 
and windows. 

The program implemented was the one-dimensional electrostatic plasma 
simulation program ESI, originally developed by Birdsall and Langdon [1985] for 
use on a Cray-1 supercomputer. While this program does not include all aspects 
of the physics that may be important for the problem we are considering, it can 
serve as a useful starting point, especially for analyzing experiments conducted 
in plasma chambers, which are less complicated than the environment of a 
tethered satellite orbiting through the ionospheric magnetoplasma. A one- 
dimensional model is clearly insufficient to deal with the combined effects of orbital 
motion and the geomagnetic field, but we should be able to use the existing 
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subroutines to develop a two-dimensional computer code. 

The original source code for ESI is written in rather obscure Fortran for a 
Cray computer. It has been translated into transparent (hopefully--at least to a C 
programmer), modular C code. The use of the C language was dictated both by 
the programmer’s preference and the greater ease with which the Macintosh 
system’s “toolbox” can be accessed with existing C implementations. The new 
version should be much easier to modify than the original Fortran version as well. 

This translation of the simulation program to a different language and a 
different system required a substantial effort, which we feel was well justified, since 
it gives us a version that we can not only use on our own computers but can share 
with other investigators using Macintosh systems. The programming effort was 
not confined to writing equivalent C code for the simulation program. It was also 
necessary to integrate this into the Macintosh environment. This required writing 
a number of i/o utility routines and graphics routines that make use of windows. 
For example, the “snapshot” plots that are periodically made and stored on 
microfiche in the original Cray version are stored as PICT files on the Macintosh 
hard disk. Software has been written to read these files and display them 
graphically in windows on the computer screen and to print them, as desired. 

We are presenting the source code for this project, as it stands. It probably 
still contains a few bugs, but it has been tested for a number of simple cases, 
including the case of counterstream tng electrons. The file-handling and plotting 
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routines are working. The simulation code has so far only been run for cases with 
periodic boundary conditions, which are not applicable to a hollow cathode model. 
The graphical display code includes a novel feature. By selecting a menu item, the 
user can make the display windows “transparent”, so that two graphs in two 
different windows (the one on top obscuring the one below it in the usual display) 
can be compared. Although we have not yet implemented it as a regular feature, 
the use of color to distinguish different particle populations has been demonstrated 
to be an easily programmed, effective analytical tool. 

To facilitate running different cases of the simulation program, while at the 
same reducing the chance of introducing errors in the input data, we have 
developed a “front end” for the program that utilizes HyperCard, the object-based 
programming environment of the Macintosh. HyperCard is easier to demonstrate 
than to describe, but the figures that follow will hopefully clarify the basic ideas. 
Figure 7 shows (with annotations) the Macintosh screen display of the first “card” 
of our plasma simulation HyperCard “stack”. All of the input parameters 
necessary to define the run, exclusive of the parameters characterizing the different 
particle species, are included on this card. That is, the number of species, number 
of time steps, etc. are displayed for the text file whose title is displayed in the 
box at the top of the card. Certain areas of the screen (card) are defined to be 
buttons (“Open”, “Save”, “Run”, “nearest grid point”, etc.). When the mouse 
button is clicked while the mouse-controlled cursor is within the defining area of 
one of these buttons, the action appropriate to that button is carried out. These 
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actions are programmable using the HyperTalk scripting language of HyperCard. 
For example, clicking on the “Run” button copies all of the input parameter values 
(as they are shown on the cards) to a selected text file (which can then be read by 
the simulation program) and launches the simulation program. The values of the 
input parameters are stored in “fields” that can be modified using the usual 
Macintosh text editing method, which relies upon selecting the text to be changed 
by means of the mouse. Only the contents of the fields can be modified in this 
way; the descriptive text is “locked”. 

The advantage of such a front-end is that the meaning of the input 
parameters can be spelled out in English, so that it is not necessary to remember 
what the equivalent variable name is in the computer program. All of the 
parameters can be displayed at once, allowing the complete picture of the computer 
run they represent to be taken in at a glance. This represents a great 
improvement over the usual alternatives, which consist of either changing the 
parameters within the source code and recompiling, making changes directly to the 
text input file (an error-prone method), or going through a tedious question and 
answer session with the computer via the screen and keyboard on every run. 
Buttons and fields are trivial to create within HyperCard, so modifications to the 
stack can be made as the program evolves with minimum programming effort. 
Figure 8 shows one of the cards that defines input parameters appropriate to a 
particular group of simulation particles. 
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size of X perturbation; 
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Figure 8 
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Details of the simulation program can be found in the Birdsall and Langdon 
reference. We now present the source code for our version, module by module. 

References 

1. Birdsall, Charles K. and A. Bruce Langdon. Plasma Physics via Computer 
Simulation , McGraw-Hill Book Company (1985). 
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I. electrostat.c 

This module includes the main control routine for the simulations. In 
“main()” the initialization routines are called and the time integration loop is gone 
through. The module also contains the event loop, which responds to mouse 
events (menu selections, etc.) whenever the program goes into the interactive 
mode. The routines to make “snapshot” and time-history graphs are included. 



include ‘plasma. h‘ 


Int my_mask = mDownMasK ♦ activMask; 

Boolean valid-event; 

FILE "output-file; 

FILE "data-read; 

double TWOPI; 

double L, dx, dt; 

double elapsed-time; 

double ael, epsi, rhoO. al, a2, eO, vsO; 

double *x. "vx, *vy; 

double rho[NGlM], phl[NGlM], e{NGIM], 

double esc{NTHl], nmslNSPMl; 

/• 

double esemlMMAXlNTHl]; 

double kes[NSPMlNTHl], pxs(NSPMINTH2]; 

V 

double *esem[MMAX]; 
double "keslNSPM], »pxs[NSPM]; 

double scratch[NGlM], fftUNGTWO]; 

double x_polnts(NGMAX]; 

double *t_polnts; 

double f_o(_v(25l; 

max_mln xv-llmitsl2]; 

double vte; 

Int mplot[MMAX]; 

Inl ng, iw; 
int nsp, ntp; 

Int It, ithl, Irho, Irhos, Iphl, le, Ixvx, Ivxvy, Ifvx; 

extern WlndowPtr draw_v/lndow{4]; 
extern PlcHandle draw_plcture[4]; 

extern Int rho-flle, rhos-flle, e_flle, phUflle, vxvy_flle, xvx_flle, 

fv_flle, hlstory_flle; 
extern long beader[I28]; 
extern Int my_ref; 

extern void do_preliminaries(); 
extern void InltlallzeO; 

extern void get-fleldsO, set_v(), bln_fv(), trace_plot(), draw_plot(); 
extern void accelerateO, advanceO, do-mouse-downO, douupdateO; 
extern void scatier_plot(), douactivateO, do_xv_plot(), do-see-updateO; 
extern void under.updateO; 

void event_loop<), recordO, finishO, show<), showlO; 
void save^fsK), save_graph(); 
double ranfO; 

Boolean RETURN; 

extern Boolean GROW, DRAG, UNDER_ONLY, INI, INF, 
extern Rect overlap_rect; 
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/. V 

mainO 

( double p-total = 0.; 
double ke_total ■0.; 
double msiNSPM], qs[NSPM). ts(NSPM]; 
double te » 0.; 
double vl, vu, vmu; 
int ins[NSPMlJ; 
register int 1th = 0; 
register int isi 
int nt; 
int nmodes; 

max-mln fv_limltst2], all_xv_lims(2j, vxvy_llmits[2]; 
do_preliminaries(ins, &nt, t>nmodes); 
tori is = 0; is < nsp; ls++) 

( initiall 2 e<lns(is}, fcinsds ^ 1], Siinslis], &qs[is], 8<ts{is}, 

&nms(ls])i 

) 

get_fields(0); 

for( is = 0; is < nsp; is+») 

( sel-vdnslls], ins(is ♦ 1], qs[ls], ms[ls], ts[ls], pxs[ls]); 

) 

SysBeepdO); 

event_loop(); 

tprlnt((output_tlle,‘ time ese P-total kel‘>; 
fprintf(output_/ile,‘ Ke2 te‘); 

whlledt < nt) 

{ iKlxvx I- 0 && (it X ixvx — 0» 

{ save-graph(xvx_file, x, vx, all_xv_lims, Ins(nsp), 

draw_window(i], &draw_plcture{l], 

•\pPhase Space ol All*, TRUE, doucv.plot); 

) 

l/dfvx 1= 0 && (it X ifvx == 0» 

{ save_fv<fv_llmlts, InsU), 13); 

) 

lf( tstO) != 0. && ivxvy !« 0 tr& (it X ivxvy •« 0)) 

{ save_graph(vxvy_file, vx, vy, vxvyJlmlts, ins[l], 

draw_wlndow(0], &draw_plcture{0], "\pVx-Vy space", 
TRUE, scatter-plot); 

) 

putotal 3 0.; 

Ke-total ° 0.; 

for( is = 0; is < nsp; ls+») 

( accelerate(ins(ls], lns[is * 1], qsUs}, ms(is], tslls], 

&pxs[lslith «li, &kes[islith]); 
pLtotal ♦= pxsdsllth ♦ 1); 
ke-total ♦« ke^isllth]; 

) 

tori is = 0; is < nsp; ls+*) 

( advanceOn^ls], ins(ls * 1], qs(ls]); 

) 

te s ke.total * eselith]; 

fprlntf(oulpuLJlle,*\nX10.3e X10,3e X10.3e X10.3c X10.3e X10.3e", 
elapsed-tlme, ese{lth}, putotal, kes(0Xith], kesUIlth], 
te); 

/•ifQth == NTH) 
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( recordO; 

W 

t-points[it] = 11* dl; 
lt++; 

elapsed-time *=■ dt; 

1th = It - Uhl; 
get-flelds(lth); 

) 

/•evenUoopO;*/ 

/*3catter_plotGc, vx, xvjlmlts, lns{l], draw_wlndow(2], 

fcdraw_picturet2], ‘NpVx versus X phase space*, TRUE);*/ 

record(nt, nmodes); 

evcnt-loopO; 

flnlshO; 

) 

/. .„«/ 

void 

cvent-loopO 

{ 

EventRecord my_event; 

WlndowPtr cvcnt-wlndow; 

Boolean valid; 

RETURN = FALSE; 
while ((RETURN) 

{ SystemTaskO; 

valid = GetNextEvent(everyEvent, &my_event); 
lf( (valid) 

{ continue; 

) 

swltch(my_event what) 

( case nullEvent: 

break; 

case mouseDown: 

lf(UNDER_ONLY I GROW) 

( break; 

) 

doi_mouse_down(8imy_event); 

break; 

case mouseUp: 
break; 

case keyDown: 
flnishO; 
case key Up: 
case autoKey: 
break; 

case updateEvt: 

event- window » (WlndowPtr)my_event. message; 
l((event-wlndow == FrontWlndowO A& 

(ietWRefCon(event_ window) < OL && INF) 

{ If(UNDER-ONLY) 

{ under_update(&my_event); 

) 

else 

{ do-see-update<8>my-event); 

} 


) 
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else 

{ douupdate<&my_«v«nt); 

UNDER_ONLY = DRAG = GROW = FALSE; 

) 

break; 
case dlskCvt: 
break; 

case activateEvt; 

do_activate(&my-event); 

break; 

case networkEvt: 
case driverEvt; 
case applEvl: 
case app2Evt: 
case appSEvt: 
case app4Evt: 
default: 

break; 

) 

) 

) 


/• V 

void 

flnlshO 

{ if(my_ref != -999) 

{ PSClose(my_ref); 

) 

ExitToShcUO; 

) 

/. V 


void 

3avc_iv(fv_Jimlts, n, nblns) 

re^stcr max_mln “fv-llmlts; 
register int n, nblns; 

{ static Int call-number = 0; 
register int save_error; 
bln_fv(vx, f_of_v, 0., 2., n, nblns); 
if(call_number 0) 

{ limit-caldx-points, f_of_v, fvJimlts, nblns); 

) 

trace-plot(x_points, f_of_v, fv_limlts, nblns, draw_wlndow(0], 

&draw_plcture{oj, *\pVeloclty Distribution*, TRUE); 
draw_plot(draw_wlndow(0], draw_plcture{0), *\pVeloclty Distribution*); 
if(call_number ~ 0) 

{ call_number*+; 

lf((save_error » wrlte(fv_flle, (char *)&header, (unslgned)512)) 

1= 512) 

{ SysBeep(20); 

) 

} 

lf((save_error = write(fv_file, (char •X*(draw_plcture{0])), 

(unslgnedX(*(draw_plcture{0]))->plcSlze))) == -1) 

{ SysBeepdO); 

) 

} 
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/. ./ 

void 

save_«raph(lplot, x, y, xy_llmits, n , graph-wlndow, graph_plc, 
graph-label, time-test, grapher) 
register Int Iplot, n; 
register double *x, *y; 
register max-min •xy-llmlts; 

WindowPtr graph-window; 

PicHandle "graph-plc; 
char ‘grap^abel; 

Boolean tlme.test; 
void (‘grapherX); 

( static int counter[14]; 
int save-crror; 

/»lt(counter[lplot] == 0 I Iplot == history-tile) HERE'S THE SCALING*/ 

{ llmlt-calctx, y, xy-limlts, n); 

) 

(•grapherXx, y, xy-limits, n, graph- window, graph-plc, graph-label, 
timeL.test); 

/*draw-plot(graph-window, *graph-pic, graph-label);*/ 
if(counter[iplot] == 0) 

{ counter[iplot}+*; 

if((save-error = wrltedplot, (char *)&header, (un3lgned)512)) 

•= 512) 

{ SysBcep(20); 

} 

) 

lt((save-crror • wrlte(iplot, (char *X**graph_pic), 

(unslgncdX(**graph_plc)->picSlze))) == -1) 

{ SysBeepdO); 

) 

) 

/. V 

void 

record(nt, nmodes) 

register int nt, nmodes; 

{ register char *modeLJabel = NewPtr(30); 
register int 1; 
max-min mode-limits[2]; 
for( 1 = 0; 1 < nmodes; i++) 

{ sprlnU(modC-label, "energy in mode Xd", 1); 

CtoPstr(mod e-label); 

save-graph(hl3tory-flle, t-polnts, esem[l), mode-limits, nt, 

draw-wlndowtl X 4), Wraw-plcture[l X 4], mode-label, 
TRUE, trace_plot); 

} 

savc-graph(hlstory-file, t-points, ese, mode-limits, nt, 

draw-wlndow(2J, 8rdraw-plcture{2], "\pelectrostatlc energy", 
TRUE, trace_plot); 

save-graph(history-tile, t_points, kes(0], mode-limits, nt, 

draw-window(3], &draw-picture{3], "\pkinetlc energy", 

TRUE, trace-plot), 

) 

/. V 

double 


ranfO 
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( double rand- value; 
double randl; 

Int first; 

randl • 1. * RandomO; 

rand-value « (randl * 32767.) / 65534.; 

return (rand. value); 

) 

/. ./ 

void 

show<x, y, n) 

register double *x, *y; 
register int n; 

{ register int i; 

ford = 0; i < n; !♦+) 

{ fprintf(output_file,*\nXe Xe', x[i], y{i]); 

} 

} 

/. ./ 

void 
showi(x, n) 

register double *x; 
register int n; 

{ register int 1; 

ford = 0; 1 < n; i++) 

{ fprintf(output_file,*\nXd ile‘,1, x[l]); 

} 

) 
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II. plasmastart.c 

This module includes routines to set up Macintosh windows and menus that 
are then available for use throughout the program. It also contains the routines 
for reading the input data file and initializing the relevant program variables. The 
first part of the actual simulation code, which creates the simulation particles 
according to the data in the input file, is also found in this module. 



♦mcluac 'plasma. n* 


WlndowRecord draw_record[4]; 

WindowPtr draw_window{4l; 

ControlHandle bars[4l2)i 
ControlHandle looKupL-controK); 

Rect drag-rect, grow_bounds; 

PlcHandle draw_plcture(4li 

Rect big_rect, biggest_rect = ( 0, 0, 756, 576); 

Rect clipper; 

int rhoJile, rh(»_fne, e^Jile, phLiile, vxvy-iile, xvx_file; 

Int fv_file, history-file; 

extern double TWOPl; 

extern double L, dx, dt; 

extern double elapsed_time; 

extern double ael, epsi, rhoO, al, a2, eO, wO; 

extern OSType my_creator, my_type; 

extern double “x, *vx, *vy; 

extern double rho(NGlM], phUNGlM], e(NGlM]; 

extern double ese(NTHl], nms(NSPM]; 

/• 

double esemlMMAXlNTHl]; 

double KeslNSPMlNTHl], pxs{NSPMlNTH2); 

*/ 

extern double *esem[MMAX]; 
extern double *keslNSPM], ‘pxstNSPMj; 

extern double scratchlNGlM], fftUNGTWO); 
extern double x_points[NGMAX], *t_polnts; 
extern double f-Of-v(25]; 
extern max_min xv_llmltst2]; 
extern double vte; 

extern int mploUMMAX]; 
extern Int ng, Iw, 
extern int nsp, ntp; 

extern int it, Uhl, irho, irhos, iphl, ie, Ixvx, Ivxvy, Ifvx; 
extern double ranfO; 

extern void accelerateO, moveJt>ars(), fill_sine-table<); 

extern FILE •output-file; 
extern FILE *data_read; 
extern Boolean RETURN; 
extern void showO; 

void fill-menusO, define_arrays(), set_defaults0, read_in_values(); 

void makc_files(), zerQ_all(), initializeO, do_even-load(); 

void doL-ordered-loadO, scrambleO, dOLj'otationO, copy_groupO; 

void add-maxwellianO, addLperturbationO, set_rho(), field-lnitO, set_v<); 

void do-inltsO, maKe^wlndowsO, set_parameters(), makc_wlndow<); 

void init-barO, swltch-controlO, do-preliminariesO, fllLxO; 
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/. 

void 

da_prellmlnarlesUns, nt, nmodes) 

Int *ins, *nt, *nmodM; 

{ Int "dummy « (lnt*X&thePort); 

dQ_init3(); 

define_arrays()i 
dummy(-63] = 331; 

TWOPI = 2. » PI; 
ntp = 100; 

It 1 Ithl = 0; 
elapsed-Ume •> 0.; 

InslO] = 0; 
zeroL^llO; 

read-in- valuesCnt, nmodes); 

maKe-filesO; 

maKe_windows(); 

L •= TWOPI; 
dx = L / ng; 
fill-x(ng ♦ 1); 
fllLsinc-UbleO; 

} 

/. ./ 

void 

dOLJnitsO 

{ InitGraK&thePort); 

InitFontsO; 

InitWindowsO; 

TElnltO; 

InltOialogs(OL); 

MaxApplZone<); 

InitMenusO; 

PlushEventst cveryEvent, 0 ); 

InitCursorO; 

MoreMastersO; 

MorcMastersO; 

MoreMastersO; 

MoreMastersO; 

MoreMastersO; 

MoreMastersO; 

StdiOLJvlaclnit(TR UE); 
sct-parametersO; 

/ill-menusO; 

} 

/. 

void 

make_windowsO 
( register int 1; 

/•blg_rect.left = 0; 
big_rect.top = 0; 

blg-rect. right » screenBlts. bounds. right; 

blgj'ect.bottom = (screenBlts.bounds.right • 10); 

blg_rcct. bottom /= 8;"/ 

big_rect = blggest-rect; 

clipper = screenBlts. bounds; 

clipper.bottom = biggest_rect. bottom; 
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/*blg_rect = screenBits. bounds;*/ 
forG = 0; 1 < 4; i**) 


{ draw_wlndow{i] = (WlndowPtr)&draw_record[l); 

make_window(draw_window{l], GRAPH_W1ND0W); 

) 

} 

/. «/ 


void 

set.parameters(> 

{ 

drag-rect = thePort->portRcct; 

SetRect<t>grow_bound3, 64, 64, thePort->portRect. right, 
thePort->portRect.bottom); 

} 

/. V 

void 

make- wlndow(ncw_ window, window-id) 
register WlndowPtr new_wlndow; 

Int wlndow_ld; 

{ ControlHandle my_scroU; 

new- window = GetNewWindow(window_id, new-wlndow, -IL); 
SctPort(new_ window); 

my_scroll = GetNewControKV-SCROLL, new_wlndow); 

SetCRefCon(my_scroll, (long)V_SCROLL); 

my-scroll = GetNewControKH-SCROLL, new- window); 

SetCRe(Con(my_scroU, (iong)H.SCROLL); 

movc_bars(ncw_wlndow); 

init-bar(new_ window, (long)H-SCROLL, 0, 50); 

init_bar(new-wlndow, (long)V_SCROLL, 0, 50); 

DrawGrowIcon(new-wlndow); 

) 

/. V 

void 

lnlt-bar<wlndow, Id, value, range) 
register WlndowPtr window; 
register Int value, range; 
long Id; 

{ ControlHandle control = lookup_control(wlndow, Id); 
il(lrange) HiliteControKcontrol, 255); 
else 

( SetCtlMlnlcontrol, 0); 

SetCtlMax(control, range); 

SetCtlValue(control, value); 

HiliteControKcontrol, 0); 

InvalRect(&<*control^>contrlRect); 

) 

) 

/. »/ 

void 

swltch-controK window, id, hillte) 

WlndowPtr window; 
long Id; 
int hillte; 

{ ControlHandle control = lookup-controKwindow, Id); 
HiliteControKcontrol, hillte); 

) 
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/. V 

ControlHandle 
lookup-controKwindow, id) 

WlndowPtr window; 
long id; 

( long label; 

ControlHandle control = «WindowPeek>window)->controlUst; , 
while (control) 

{ label = GetCRefCon(control); 
if (label == id) break; 
control = (•control)->ncxtControl; 

) 

return control; 

} 

/. ./ 

void 

fill_inenus() 

{ char *apple_name = NewPtr(20); 

MenuHandle menu; 

‘applcLJiame = (char)l; 

‘(apple^jiame « 1) = (char)appleMark; 
menu = NewMenu(APPLE_MENU, apple_name); 

AppendMenu(menu, *\pAbout Plasma ...''); 

AppendMenu(menu, *\p(-*); 

AddResMenu(menu, DRVR); 

InsertMenu(menu, 0); 

In3ertMenu(GetMenu(PILE-MENU), 0); 
InsertMenu(QetMenu(EDIT-MENU), 0); 
lnsertMenu(menu - (>etMenu(PLASMA-MENU), 0); 

CtoPstr(apple_name = •) Shrink to Pit"); 
apple_name(2] = (char)checkMark; 

AppendMenu(menu, apple-name); 

OrawMenuBarO; 

Dlsableltem(GetMenu(£01T_MENU), 0); 

) 

/. ./ 

void 

deflne_arrays() 

( register int 1; 

X = (double •)NewPtr(slzeof(double{NPAR])); 
vx = (double •)NewPtr(sl 2 eof(double(NPAR])); 
vy = (double •)NewPtr(sizeof(double{NPAR])); 
for( 1 » 0; 1< NSPM; !♦♦) 

{ pxs(i] ° (double *)NewPtr(NTH2 * sizeof(double)); 
kes[i] > (double *)NewPtr(sizeof(double{NTHl])); 

) 

t_points = (double *)NewPtr(sizeof(double{NTHl])); 
for( 1 = 0; 1 < MMAX; !♦♦) 


{ esemCi] = (double •)NewPtr(NTHl • slzeof(double)); 

} 

} 

/. ./ 


void 

set_defaults(nt) 
int *nt; 
register int i; 


{ 
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nsp = 1; 

L = TWOPIj 
dt = .2; 

•nt = 10; 
epsl = 1.; 
ng = 32; 
iw = 2; 

al = a2 = eO s wO = 0.; 

Irho = irhos = tphl = Ic = ixvx = Ivxvy = ifvx = 20; 
ford = 0; i < MMAX; !♦♦) 

I mploUil = 0 ; 

) 

) 

f* ./ 

void 

read_ln_values(nt, nmodes) 

register int *nt, ‘nmodes; 

( register int i; 

char ch_read =(charX>; 

data_read • fopen(*simulation data*, *r'); 

while(ch-read 1= (*\n*XOl) 

( fscanf(data_read,*Xc*,&cb_read); 

) 

fscanf<data.read,*Xd\nXd\nXlf\nXlf\nXd\n*, Ansp, nt, Mt, AL, gng); 
fscanf(data_read,*Xlf\nXlAnXd\nXd\nXd\nXd\nXd\n!Cd\nXd\nXd\n*, 
&epsi, irhoO, &iw, &lrho, Airhos, &iphi, &ie, &ixvx, Wvxvy, 
Aifvr); 

fscanf(data_read, *Xlf\nXlAnXlAnXlAnXd\n*, &al, &a2, 4«0, &wO, 
nmodes); 

ford = 0; i < ‘nmodes; i++) 

{ rnploUi] • i; 

} 

) 

/. ./ 

void 

make_files() 

{ register int open_mode = 0_RDWR + 0_CREAT + 0-BlNARY; 
output_flle = fopen(’run_flle", *w+“); 
my_type = 'PLOT'; 
my_creator = 'PLAS'; 
ifdrhos 1= 0) 

{ rhos_fiie = open(*smoothed rho.plot", open_mode); 

) 

ifdrho I" 0) 

{ rho-file = open<*rho.plot‘, open_mode); 

} 

ifde (• 0) 

{ e_file = openCelectrlc field. plot’, open_mode); 

) 

Ifdphi 1= 0) 

{ phi-file = openCphi.plof, open_mode); 

) 

Ifdvxvy 1= 0) 

{ vxvy_flle = open(*vxvy.plof, open-mode); 

) 

ifdxvx 1= 0) 
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( xvx_fUc = openCxvx.plof, open_mode); 

) 

ifdfvx != 0) 

{ - open(*f(v).p!ot*, open-mode); 

) 

hlstory_/lle = opcnChlstory.plof, open_mode); 

) 

/. 

void 

zero-allO 

{ register Int 1, J; 

for( 1 = 0; 1 < NSPM; !♦♦) 

( pxstilO] = 0.; 

for( J = 0; J < NTHl; J++) 

( KestiU) = 0.; 

pxsllU + 11= 0.; 

> 

} 

for( 1 = 0; 1 < NTHl; 1++) 

( esell] = 0.; 

1 

for( 1 = 0; 1 < MMAX; 1++) 

{ mploUi] = 0; 

forO = 0; J < NTHl; J++) 

{ esetndUl = 0.; 

} 

} 

) 

/. 

void 

fill_x(n) 

register Int n; 

{ register Int 1; 

ford = 0; 1 < n; 1++) 

{ x_points[l] = 1; 

} 

} 


•V 


/. V 

void 

Inltlallzedll, 112, m, q, t, nm) 
double *m, *nm, *q, *t; 
register Int ill, *112; 

( double Lg, ddx, wp, wc, vtl, vt2, vO, xl, vl, thetax, thetav, qm; 
int ngr, n, nig, nv2, mode; 

Boolean sKip = FALSE; 

fscanf(data_read,‘']|d\nXlf\nXlAnXlf\nXd\n*; An, Awp, &wc, Aqm, Anlg); 
fscanf(data_read,*XlAnXlAnXlAn!M\n", &v6, ivtl, Avl2, Anv2); 
fscanf(data_read,'Xd\nXlAnXlAnXlAnXlAn’, Amode, Axl, Avl, 

Athetax, Athetav); 

Ifdll == 0) 

{ vte = vtl ♦ vt2; 

) 

*t = tan( -wc • dt • .5); 

*112 = ill + n; 

•q = L * wp * wp / (n * qm); 
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•m = *q / qm; 

•nm = n • (*m); 

ngr = n / nig; 

Lg = L / nig; 
ddx = L / n; 

da-even-load<lll, ngr, vO, ddx); 

«(ORDERED) 

{ do-ordered_load(ill, *112, n, nv2, ngr, vt2, &sklp); 
If(lsklp) 

{ scramble(Lg, ddx. 111, ngr); 

) 

) 

ifCMAGNETIZEO && (skip) 

{ dOLj'otalionlngr, 111, Lg); 

} 

IfCANOTHER) 

{ copy_group(m, ngr, n, Lg, wc); 

) 

If(RANDOM) 

( add_maxwellian<n. 111, vtl, wc); 

} 

add_pcrturbatlon(n, 111, mode, xl, vl, thetax, thetav); 
sct_rho(ill, *112, *q, n * (»q) / L); , 

/*3howl(rho, ng ♦ 1);*/ 

) 

/. 

void 

do^evea-loaddll, ngr, vO, ddx) 
register Int ill, ngr; 
double vO, ddx; 

( register int i, il; 
double xO; 

tor( 1 = 0; 1 < ngr; !♦♦) 

{ 11 = I + ill; 

xO = (1 » 1. ♦ .6) » ddx; 
x[il] = xO; 
vx[ll] = vO; 

} 

) 

/. 

void 

da_ordered_load(ill, 112, n, nv2, ngr, vt2, skip) 
int ill, 112, nv2; 
register int n, ngr; 
double vt2; 

Boolean *sklp; 

( register Int 1, 11, J; 

double vmax, dv, w, wnv2, tv, df; 

vmax = .5 * vt2; 

dv = 2. * vmax / (n - 1); 

wnv2 = 1.; 

xlm] - 0.; 

ford = 1; 1 < n; 1++) 

{ w » (d - .5) * dv - vmax) / vt2; 
lf(nv2 1= 0) 

( wnv2 = pow<w, nv2); 


•*/ 





) 

fv = wnv2 • exp( -.5 • w • w); 

11 = i ♦ HI; 

x(ll] » xfll - 1] ♦ my_max(fv, .0); 

} 

<« = x[ll] / ngr; 

11 = 111 ; 

J = 111; 

for( 1 « 0; 1 < ngr; 11++) 

{ fv • (1 ♦ .5) • df; 

whlledv >= x[J + 1]) 

( J++; 

lf(J > 112 - 2) 

{ ‘skip = TRUE; 
return; 

) 

) 

w = dv * (J - til ♦ (fv - x[J]) / (x[J ♦ 1] - xlj])) - vmax; 
vx(ll] ♦= w; 

) 

) 

/• 

void 

3cramble(Lg, ddx, ill, ngr) 
double Lg, ddx; 
register Int ill, ngr; 

{ register int i, ii; 
double xs <• 0.; 
double xsi; 

for( 1 » 0; 1 < ngr; !♦♦) 

{ il = 1 ♦ ill; 

xCll] = xs ■ Lg ♦ .5 • ddx; 

/» wrlteO; »/ 
xsi >= 1.; 
while(xs >= 0.) 

{ xsi «= .6; 
xs -= xsi; 

) 

xs +- 2. • xsi; 

) 

) 

/. 

void 

da_rotation(ngr, ill, Lg) 

register int ngr, ill; 
double Lg; 

( register int i, il; 
double w, theta; 
for( 1 = 0; i < ngr; !♦+) 

{ il = ill ♦ i; 
w = vxlilj; 

theta = TWOPI ‘ x(ll] / Lg; 
vxlll] = w ■ cosUheta); 
vytll] = w • sln(theta); 

) 
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/. ./ 

void 

copy_group(ill, ngr, n, Lg, wc) 
int ill, ngri 
register int n; 
double Lg, wc; 

{ register int i, J, il, 12, 
double xs s 0.; 
for (i = ngr; i < n; l ♦= ngr) 

( xs ♦= Lg; 

for( J = 0; J < ngr; >♦> 

{ il = J ♦ ill; 

12 = 11 ♦ i; 
x(i2] = x(ll] ♦ xs; 
vx[i2] = vxtil]; 
ifCMAGNETIZED) 

{ vyfl2] = vy(ll]; 

} 

) 

) 

} 

/. V 

void 

sdd-maxwellian(n, ill, vti, wc) 
register int n, ill; 
double vtl, wc; 

{ register Int 1, 11, J; 
for( 1 “ 0; 1 < n; !♦♦) 

{ 11 = ill ♦ 1; 

for( J = 0; J < 12; J++) 

{ If(MAGNETlZED) 

( vy(ll) ♦= vtl • (ranfO -.5); 

) 

vx[ll] += vtl ■ (ranfO -.5); 

) 

) 

} 

/. ,-«/ 

void 

add_perturbatlon(n, ill, mode, xl, vl, tlietax, thetav) 
register int n, ill, mode; 
double xl, vl, thetax, thetav; 

{ register int i, 11; 
double theta; 
for( 1 = 0; 1 < n; !♦♦) 

( 11 = ill ♦ 1; 

theta ■ TWOPl * mode * xtil) / L; 
xtll) xl * cositheta ♦ thetax); 
vxtllj ♦= vl • sinCtheta ♦thetav); 

) 

) 

/. ./ 

void 

seLj'hoOl, iu, q, rhos) 
register int 11, iu; 
double q, rhos; 



( register int 1, J; 

double dxi, xn, drho, qdx; 
qdx = q / dx; 
dxi • 1. / dx; 
xn = ngi 
lf(FIRST_GROUP) 

{ for( 1 = 0; i < ng; !♦♦) 

{ rholi] = rhoO; 

> 

rhoing] = 0.; 

) 

/• show<x-points, rho, ng ♦ I);V 
rhoO -= rhosi 
tor( 1 = 0; 1 < ng; 1++) 

( rho(l] -= rhos; 

} 

/* show<x-polnts, rho, ng ♦ l);*/ 
swltchGw) 

( case ZERO.OSOER: 

for( i = il; i < lu; i**) 

{ x[l]*=dxl; 

»(xUl < 0.) 

{ x[ll ♦= xn; 

} 

lf(x[l] >= xn> 

{ x[i] -= xn; 

} 

j « x(i] ♦ .5; 
rho(J} ♦= qdx; 

} 

break; 

case MOMENTUM: 
case ENERGY: 

for< i = U; i < lu; !♦♦) 

{ x[i) •= dxi; 

l/(xll) < 0.) 

( x[ij ♦= xn; 

) 

if(x(il >= xn) 

{ x(l] -= xn; 

) 

j = xm 

drho = qdx • ( x[l] - J); 
rho(J] ♦= (qdx - drho); 
rho(J ♦ 1] ♦= drho; 

} 

break; 

) 

/• show(x_points, rho, ng + !);•/ 

) 

/. 

void 

flelcLJnit(sni, ksql, ng2) 

register double *sm, *ksqi; 
register int ng2; 



I register int 4 
double kdx3j 
for( 1 = 0; 1 < ng2; !♦+) 

{ kdx2 = (PI / ng) ■ (1 ♦ 1); 

sm[i] 3 exp(al * pow(sin(Kdx2), 2.) - a2 * pow(Un(kdx2), 4.)); 
KsqUi] = (pow<9mtll • dx / (2. • sln(Kdx2)), 2.) > / cpsl; 

) 

} 

/. 

void 

seL.v(il, iu, q, m, t, p) 

double q, m, t, *p; 
register Inl 11, lu; 

{ double dtdx ° dt / dx; 
double c, s, vxx; 
register Int 1; 
ll(NEED_ROTAT10N) 

{ c = 1. / sqrtd. ♦ t • t); 
s = c • t; 

for ( 1 • 11; 1 < iu; !♦♦) 

{ vxx = vx[l]; 

vxtl] = c • vxx ♦ s • vytD; 
vy[l] = -3 ■ vxx ♦ c • vytO; 
v^ll •= dtdx; 

) 

} 

for < 1 = 11; 1 < iu; 1++) 

{ vxfl] •= dtdx; 

) 

acceleratedl, lu, -.5 • q, m, 0., p, scratch); 

} 
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III. plasmafiles.c 

This module consists of the routines to create, read, display, and print 
graphics files: snapshots made at selected intervals or time histories made the end 
of the simulations. 
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^Include 'pldsma.tr 
int dirty; 

Size file-size; 
long header(126]; 

Point sfg_where = {90, 82); 

SPReply reply; 

Str255 picture_name; 
int my_ref = -999; 

PlcHandle file_picture; 

extern PILE *output_file; 
extern Boolean MOKE_PLOTS; 
extern Rect big_rect, clipper; 

extern void show.scrollsO, hide_3crolls(), erase_grow{), dOL.wlndow_care<); 

void do_new(), dOL^ve-asO, dOL.open(), do_save(), show_error(); 
void doureadO, douprint_close<), handle_error(); 

/• - - •/ 

do_print() 

{ GrafPtr savei^af, save^print; 

WindowPtr the^-window = ProntWlndowO; 

TPPrPort print-port; 

THPrint my_print; 

TPrStatus pL.status; 

TSetRslBlk lmage_set; 

Boolean valid; 
int ‘prlnt-rect; 
int i; 

PicHandle draw-picture = (PicHandle)my_abs(GctWRefCon(ProntWindow())); 

GetPort(«rsave_graf); 

image_set.iOpCode = setRslOp; 

image_set.iXRsl - 144; 

image_set.lYRsl = 144; 

PrOpenO; 

IfCPrlntErr 1= noErr) 

{ return PrintErr; 

) 

my_prlnt » (THPrint)NewHandle<5lzeof(TPrlnt)); 

PrlntOefault(my_prlnt); 
valid = PrValidate(my_print); 
imagc-set.hPrint = my_print; 

PrGeneral(&image_set); 
ifdmage-set.iError 1= 0) 

( SysBeepdO); 

SysBeep(10); 

) 

print-port = PrOpenDoc<my_print, OL, OL); 
ifCPrintErr 1= noErr) 

{ return PrintErr; 

} 

GetPort(&saveL-prlnt); 

SetPort(the_ window); 
hldeL-ScrollsUheL-wlndow); 
eras€L 4 row<thc_ window); 

SetPort(save_prlnt); 

valid = PrJobDialog(my_prlnt); 
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ifOvalid) 

( SatPort(tha_wlndow); 

SysBeepdO); 

show_scrolls(the_ window); 

DrawGrowlcon(the_wlndow>; 

SetPort(save_graf); 

return; 

) 

valid = PrStlDialog(my_print); 
ifdvalid) 

{ SctPort(the_window); 

show_3croll3(the_window); 

DrawGrowIcon(the_window); 

SetPort(save_graf); 

return; 

) 

PrOpcnPage(prlnt_port, OL /“iblg-rect*/); /* this choice means no scaling •/ 
IfCPrintErr 1= noErr) 

{ do_print_closc(prlnt_port, save_graf); 

return PrintErr; 

) 

print_rcct = (Int *)4(prlnt_port->gPort.portRect); 
for( i = 0; 1 < 4; i++) 

{ fprintf(output_file, *\nlM*, *(print_recU+)); 

} 

DrawPicture(draw_picture, fc(prtnt_port->gPort.portRect) 

/•&big_rect*/); 

if(PrintErr !• noErr) 

( dQ_print_closc(prlnt_port, save_graf); 

return PrintErr; 

) 

PrClosePage<prlnt_port); 

PrCloseDoc(prlnt_port); 

if«*my_print)->prJob.bJDocLoop == bSpoolLoop) 

{ if(PrlntErr l« noErr) 

{ PrCloseO; 

SetPort(save-graf); 
return PrintErr; 

) 

PrPicFlle(my_print, OL, OL, OL, ipustatus), 

) 

PrCloseO; 

SetPort(the_ Window); 
show_3crolls(the_ window); 

DrawGrowlcon(the_ window); 

SetPort(save_graf); 

) 

/. «/ 

void 

douprlnt_close(prlnt_port, the_save) 

TPPrPort prlnt_port; 

GrafPtr the^save; 

{ PrClosePage(prlnt_port); 

PrCloscDoc<print_port); 

PrCloseO; 

SctPort(the-save); 
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) 

/. ./ 

void 

handle_error(the_error) 

Int the_error; 

( 

) 

/. V 

void 

do_new<) 

{ /•SetWTltle<draw_window(0],"\pUntiU«l"); 

ShowWindow(draw_window{0]); 
dirty = 0;V 

} 

/. V 

void 

do_savc_as() 

{ register WindowPtr plot-window = FrontWlndowO; 

register int i; 

Str255 theL.volume; 
long file-count; 

register PlcHandle plot-picture; 
long header-X^ount = 512; 
if(plot- window == OL) return; 

plot-picture = (PlcHandle)GetWRefCon(plot- window); 
file-size = file-count = GetHandleSizefplot-picture); 
plcturc-name{0] = 0; 
hlde-Scrolls(plot_ window); 
erase-grow(plot_ window); 

SFPutPlle<sfg_where, “\pSave File As plcture^name, OL, Areply); 
3how_scrolls(plot_window); 

DrawGrowIcon(plot-window); 
if (reply good) 

{ lf((l = FSOpenfreply.fName, reply. vRefNum, &my-ref» 1= fnfErr) 

{ ; 

) 

else 

{ 1 = Createfreply.fName, reply. vRefNum, '7?77‘, 'PICT'); 

i = FSOpen(reply.fName, reply. vRefNum, &my_ref), 

) 

SetWTltlefplot- window, reply. fName); 

i = FSWrlte(my-ref, &header-count, ‘header); 

HLocK(ploL-picture); 

1 = FSWrite^y-ref, &fiieL.count, ‘plot-picture); 

HUnlock(plot-picture); 

lf(l 1= noErr) 

{ show-error(l, 30); 

} 

else 

{ lf((i = GetVol((StringPtr)theL-volume, Areply. vRefNum)) 

1= noErr) 

{ show-errorO, 40); 

) 

lf((i = FSClose(my-ref)) != noErr) 

{ show-error(i, 50); 

} 
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lf((l = PlushVol((StrlngPtr)theL.volume, reply .vRefNum)) 
Is noErr) 

{ show_error<l, 60 ); 

) 

) 

) 

) 

/. V 


void 

show_error(typeL.error, which_call) 
int type^crror, which^all; 

( SysBeep(which_call * 10); 

fprintf(output_file, *\nerror ° Ed on call number Ed", type-error, 
which-call); 


) 

/. 

void 

douopen(new_ window, new-picture) 
WlndowPtr new_window; 
PicHandle “new_picture; 

( SPTypeUst my_types; 

WindowPtr first- window; 
register int i; 
long header-count = 512; 
MORE-PLOTS = FALSE; 


my_types(0] = 'PLOT'; 
if(my_ref != 999) 

( 1 - FSClose(my_ref); 

my_ref = -999; 

} 

IfWfirst- window = ProntWindowO) 1= OL) 

{ hidc-scrolls(first_wlndow); 

erase-grow(first-window); 

) 

SFGetFile(sfg-where, *\p*, OL, 1, my_types, OL, fcreply); 
if(first_window 1= OL) 

{ show_scrolls(first_window); 

Draw<5rowlcon<flrst_window); 

) 


»/ 


if (reply. good) 

{ if((l = FSOpen(reply.fName, reply. vRefNum, &my-ref)) 1= noErr) 
{ show-error(l, i); 
return; 


} 

if((i = FSRead(my_ref, &header_count, Aheader)) 1= noErr) 
{ show_error(l, 2); 
return; 


) 


doLj'ead(new_window, new_picture, reply, my_ref); 

} 


else 

{ do-wlndow-careO; 

} 


} 

/• 


void 


-V 
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do_read(new_wlndow, new.picture, thCLj'eply, the_ref) 
register WindowPtr new.window; 

PicHandle *new_picture; 

SFReply thcLj'eplyi 
register int the_refi 
( register int i; 
int J; 

long piCLXOunt « 2L, 

WindowPtr first-window; 

MORE-PLOTS = TRUE; 

if((i c PSRead(theLj'ef, &piCLXOunt, &J)> 1= noErr) 

{ if(i = eofErr) 

{ if((first-Window = ProntWindowO) 1= OL) 

{ hideL^rolls(first-wlndow); 

erasei^owffirst- window); 

} 

NoteAlertfENOALERT, OL); 
if(first-window 1= OL) 

{ show-scrolls(flrst-wlndow); 

Draw6rowlcon(first_window); 

} 

1 = PSClose(the-ref); 
my-ref = -999; 

MOREJ>LOTS = FALSE; 

) 

return; 

} 

•new-picture - (PtcHandle)NewHandle((long)J); 

HLock(*new-plcturc); 

if((l = SetFPos(the-ref, fsFromMark, -2L» 1= noErr) 

{ show_error(i, 4); 

return; 

) 

pic-count = Uong)J; 

lf((i = FSReadCthcLj'ef, &plc-count, ••new_plcturc)) 1= noErr) 

{ if(i • eofErr) 

{ NoteAlerUENDALERT, OL); 
i = FSCloseftheLJ-ef); 
my_ref « -999; 

MORE-PLOTS = FALSE; 

} 

} 

SetWRefConOiew-wlndow, Oong)*new_picture); 

SetPortOiew- window); 

ClipRectf&clipper); 

EraseRectf&clipper); 

SetWTitle(new-window, &the_reply.fName); 

SelectWindow(new_wlndow); 

dOL.wlndow_care<); 

ShowWindow(new-window); 

DrawPlcture(*new_plcture, &big_rect/*&new_window->portRect*/ ); 
HUnlock(*new-picture); 

/“DrawControls(new-wlndow); 

DrawGrowIcon(new-wlndow);*/ 

/■ THE FOLLOWING CODE IS BASICALLY RIGHT FOR DEFINING A BITMAP 
PICTURE OF THE GRAPH. BUT IT ISN'T INTEGRATED WITH THE REST 
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OP THE PROGRAM. LOTS OP BIZARRE DRAWING RESULTS IN TRANSPARENT 
MODE. 

•new_plcturc = 6penPicture(&big-rect); 

CopyBUs(&new_window->portBlts, &new_window->portBits, &big_rect, 
&big_rect, srcOr, OL); 

ClosePictureO; < 

SetWRefCon(new_window, Oong)*new_picture)i 

V 

) 

void 

dcusaveO 

{ 

) 
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IV. plotting.c 

This module contains routines to make scatter and line plots, which are 
stored as “pictures” available for display through the routines found in 
plasmafiles.c. 
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^include *plasma.b‘ 

extern Reel big_rect, clipper; 
extern PILE ‘output-file; 
extern double elapsed-tlme; 

extern double L, dx, dt; 

void scatter-plotO, llmlt-calcO, trace_plotO, llneUtO, move_lt(); 
void map-pointO, make_pointO, bin_f^), draw_plot(), dOLjcv_plot(); 
void xv_limit_calc(), dQ_wlndow_care<>; 

/. ./ 

void 

scatter-plotCx, vx, xv_limlts, n, plot-window, plot-picture, plot-label, 
time-test) 
register int n; 
register double ‘x, *vx; 
register max-min ‘xv-limits; 

WindowPtr plot-window; 

PicHandle ‘plot-picture; 
char ‘plot-label; 

Boolean time-test; 

( register int i; 

Rect plot-rect; 

double deltax, deltay; 

char ‘e-string = NewPtr(256); 

SelectWlndow(plot_ window); 

ShowWlndow(plot_ window); 

SetPort(plot-Window); 

ClipRecU&clipper); 

EraseRect(&clipper); 

plot.rect.top < 11; 

plot-rect.bottom ° 261; 

plot-rect. left = 31; 

plot-rect. right = 481; 

limlt-calc<x, vx, xv-limits, n); 

deltax s xvjlmits(0].fmax - xv-limits[0].fmin; 

deltay - xv-llmlts[l}.fmax - xv_limlts(l].fmin; 

‘plot-picture - OpenPicture<&big-rect); 

ShowPenO; 

PrameRecU&plot-rect); 
ford = 0; 1 < n; !♦♦) 

( make-polnt((x[i] - xv_limits[0].fmin),(vxll] - xv-limlts[l].fmin), 

deltax, deltay); 

) 

MoveTo(180,295); 

DrawString(plot-label); 

If(time-test) 

{ sprintf(e-strlng, 'time = X.21f, elapsed-tlme); 

CtoP8tr(e-strlng); 

MoveTo(350, 276); 

DrawString(e^string); 

) 

ClosePictureO; 

DisposPtr(eL^trlng); 

SetWRefCon(plot- window, (long)‘plot-picture); 



dQ_wlndow_care(>; 

DrawControls(plot_window)i 

Draw<>rowlcon(plot_wlndow); 

) 

/. 

void 

do-xv_plot(x, vx, xv-limits, n, plot-window, plot-picture, plot-label, 
time^test) 
register Int n; 
register double *x, *vx; 
register max_Riin *xv_limlts; 

WindowPtr plot-window; 

PicHandle *plot_picture; 
char *plot-label; 

Boolean time-test; 

( register int 1; 

Rect plot-rect; 
double deltax, deltay; 
double dxdt = dx / dt; 
char “e-string • NewPtr(256); 

SelectWindow<plot- window); 

ShowWindow<plot- window); 

SetPort(plot- window); 

ClipRect(&cllpper); 

EraseRect(&clipper); 

plot-rect. top = 11; 

plot-rect.bottom = 261; 

plot-rect. left > 31; 

plot-rect. right = 481; 

xv-limit_calc<x, vx, xv-limtts, n); 

deltax ° xv_limits(0].fmax - xv-limitstO].fmin; 

deltay = xv-llmlts(l].fmax - xv-llmltstl).fmin; 

•plot-picture = OpenPlcture(&blg-rect); 

ShowPenO; 

FrameRect(&plot-rect); 

ForeColor<redColor); 
ford = 0; 1 < n; 1++) 

{ 1/(1 >127) 

{ ForeColor(blueColor); 

) 

maKe-point((x[l] - .5 • vx[l) • dx - xv-llmltslO).fmln), 
(vx[i] • dxdt - xv-limits(l].fmin), 
deltax, deltay); 

} 

ForeColor(blackColor); 

sprlntf(e-string, “y min = *.2le", xv_llmlts(l].fmln); 
CtoPstr(e-strlng); 

MoveTo(38, 276); 

OrawString(e-string); 

sprintfCcL^trlng, "y max = X.21e‘, xv-llmltstl).fmax); 
CtoPstr(e-string); 

MoveTo(200, 276); 

DrawStrlng(eLjtring); 

MoveTo(180,295); 

DrawString(plot-label); 

if(time^test) 



( spnnu<e_strmg, *ume « eiapsea-ume)^ 

CtoPstHft-string); 

MovcTo(350. 276); 

DrawSlring(e_9lring); 

) 

ClosePlctureO; 

DisposPtr(e_string); 

SetWRefCon(plot_window, (long)*plot-plcture); 
dOL.window_jcare<); 

DrawControls(plot- window); 
Draw6rowlcon(ploL.wlndow); 

} 

/* 

void 

draw_plot(plot_wlndow, plot_plcture, plot-label) 
register WlndowPtr plot-wlndow; 

PlcHandle plot-plcture; 
char *plot_label; 

{ /•SelPort(plot-wlndow); 

CllpRectC&clipper); 

ErascRecU&cllppcr);*/ 

SetWRefCon(plot_window, (long)plot_plcture); 
SetWTltle(plot_window, plot-label); 

} 

/. 

void 

limlt-calc(fl, 12, f_llmlts, n) 
register double *11, *12; 
register max_mln “f-llmlts; 

( register int 1; 

double 111 s l.elO; 
double flu = -l.elO; 
double 121 » l.elO; 
double 12u » -l.elO; 
for ( 1 = 0; 1 < n; 1++) 

( if(fl[i] < 111) 

{ fll = fill]; 

) 

iKfUl] > flu) 

{ flu = fl[l]; 

) 

if(f2[i] < f21) 

{ f21 = f2[l]; 

} 

lf(f2[l] > f2u) 

{ f2u = f2[lj; 

) 

) 

f-limlts(0].fmax ° flu; 
f-limits[0].lmin = fll; 
f-limits[l].lmax = f2u; 
f_llmltsll].fmln = f21; 

} 

/» 

void 

xvJlmlt-calcffi, f2, f-llmits, n) 



register double *fl, ■f2; 
register msx-mln *f_llmlts; 

{ register Int 1; 

double fU - l.elO; 
double flu = -l.elO; 
double f21 c l.elO; 
double f2u > -l.elO; 
double gl, g2; 
double dxdt dx / dt; 
for ( 1 = 0; 1 < n; 1++) 

{ gl = fl[l] - .5 « f2[l] • dx; 
g2 = f2[ll » dxdt; 
lf(gl < fll) 

{ fll = gl; 

) 

lf(gl > flu) 

( flu = gl; 

} 

«(g2 < 121) 

{ f21 - g2; 

) 

lf(g2 > f2u) 

{ f2u = g2; 

) 

) 

fJlmltrfO],fmax = flu; 
f_llmltsIO].fmln * fll; 
fjllnilts[l].fmax - f2u; 
fJllmlts(l].fmln = f21; 

) 

/. 

void 

trace-ploUx, y, xy_llmlts, n, plot-wlndow, plot_plcture, plotJabel, 
time-test) 
register double *x, *y; 
register max_mln *xy-llmlts; 
char *plot-label; 
register Int n; 

WindowPtr plot-wlndow; 

PicHandle ■plot_plcture; 

Boolean tlme.test; 

{ register Int 1; 

Rect plot_rect; 

double deltax, deltay; 

char *e_string = NewPtr(256); 

SelectWindow(plot_wlndow); 

ShowWindow(plot_ window); 

SetPort(plot_wlndow); 

ClipRecti&clipper); 

EraseRect(t«lipper); 
plot-rect.top = 11; 
plot-rect. bottom - 261; 
plot_rect.left = 31; 
plot-rect. right = 481; 

deltax - xy-Umlts[oi.fmax - xy-limltslO].fmin; 
deltay = xy-limits(l].fmax - xy-limits[l].fmin; 
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•ploUplcturc = OpenPlcturc(*<Wg_rect); 

ShowPenO; 

FrameRecU&plot-rect); 
move_it(x[0] - ry_limUs(0].fmin, 

y(0] - xy_limitsn].fmln, delUx, delUy); 
ford = 0; 1 < n; !♦♦) 

{ UneLJt{x[il - xy-limits[0].fmin, 

^i] - xy_limits(l].fmin, deltax, deltay); 

} 

/*move_lt(x[6] - xy-limits(0].ftnin, 0., deltax, deltay); 
lineJt(x[6] - xy_limlts[0}.fmln. deltay, deltax, deltay);*/ 
sprintf(e_string, *y min ° X.21e‘, xy-limlts[l].fmln)i 
CtoPstKe-string); 

MovcTo<38, 276); 

DrawString(e-String); 

sprintf(e_string, *y max = X.21e", xy_llmlt5[l].fmax); 
CtoPatr(e_strlng); 

MoveTo(200, 276); 

DrawString(e_string); 

MoveTo(180,295); 

DrawStrlng(plot_label); 

If(tlme-test) 

( sprintf(eL^tring, "time = X.21f‘, elapsed-time); 
CtoPstrifiL-strlng); 

MoveTo(350, 276); 

DrawString(c-string); 

) 

ClosePlctureO; 

Diapo3Ptr(e_5trlng); 

SetWRefCon(plot- window, (long)*plot_picture); 
dOL.window_care(); 

DrawControla(plot_ window); 

DrawGrowIcon(plot_window); 

) 


/. «/ 

void 

move_it(x, y, dx, dy) 

double X, y, dx, dy; 

( Int h, v; 

mapi-poinUx, y, dx, dy, &h, &v); 

MoveTo(h, v); 

) 

/. »/ 

void 

line_lt(x, y, dx, dy) 

double X, y, dx, dy; 

{ int h, v; 

map-point(x, y, dx, dy, &h, dv); 

LlneToOi, v); 

} 

/» */ 


void 

map_polnt(x, y, deltax, deltay, h, v) 
double X, y; 
double deltax, deltay; 
reglater int ‘h, *v; 
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( double tst; 

U(deltax > l.e-46) 

{ X /» deltax; 

■h = ( X * 450); 

) 

else 

{ *h = 225; 

) 

•h ♦= 31; 

iKdeltay > l.e-45) 

{ y /= deltay; 

*v = 250 - (y » 250); 

) 

else 

{ »v = 125; 

} 

»v ♦= 11; 

) 

/. ./ 

void 

make_pointG(, y, deltax, deltay) 
double X, y; 
double deltax, deltay; 

( Rect pt_rect; 

Int h, v; 

register Int r = 1; 

map-poinUx, y, deltax, deltay, &h, &v); 

pt_rect.top • V - r; 

pLrect.left = h - r; 

pt_rect. bottom = v ♦ r; 

pt_rect. right = h ♦ r; 

PalntOval(&pt_rect); 

} 

/. V 

void 

bin_fv(v, f.dist, vcenter, dv, n, nbins) 
double vcenter, dv; 
register double *f_dlst, *v; 
register int n, nbins; 

{ register int 1, J; 

double space = 2. ■ dv / nbins; 
double compare, coml; 
coml = vcenter - space • (.5 • nbins -1); 
for(J = 0; J < nbins; J++) 

{ f_dl3t[J] = 0.; 

) 

for { 1 = 0; 1 < n; 1++) 

{ compare = coml; 

for(J = 0; J < nbins; J*+) 

{ if(v[l] <= compare) 

{ f_distU) ♦= 1.; 
break; 

} 

compare ♦= space; 

) 

) 


) 
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V. response.c 

This module includes the routines that take care of all interactive mode 
events. In this mode menu selections (e.g. examining plot files) are executed and 
graphics viewing windows are moved and selected. The code for making a window 
“transparent” is found in these routines. 



Page 94 


*inciuae ’piasma.n* 

•daflns ENOUGH 3 
•define NOT-ENOUGH 2 

int DESKOPEN = 0; 

extern WlndowPtr draw_wlndow(4]; 
extern PicHandle draw_plcture[4]; 
extern PILE *output_flle; 

extern int my_ref; 
extern SFReply reply; 

extern Reel drag_rect, growjaounde, blg_rect; 
extern Boolean RETURN; 

extern void dcunewO, do-save^sO, doLU>penO, do-saveO; 
extern void do_readO, handleL_errorO, swltch-controlO; 
extern dO-prlntO; 

void dO-controlsO, dounenuO, dcuaboutjnenuO, douJllesO; 
void do_mou9e_down0, douupdateO, dcuexamlneO, move-barsO; 
void grow_wlndow(>, invaLgrowO, do_activateO, show_scrollsO; 
void hlde_scrolls(), dOL.wlndow_care(), erase_gro^), valld-growO; 
void toggle^seeO, combineL.pictures(), valid-scrollsO, grow.updateO; 
void valid-rectO, defineL4P'0w<), dO-see-updateO, under.updateO; 

Boolean MORE_PLOTS, INI, INF, see_through(), GROW, DRAG; 

Boolean UNDER_ONLY; 

PicHandle combo-pic; 

Rect savej"ect, overlap_rect, save^row, save_v, save_h; 

WlndowPtr top_window, under.window; 

RgnHandle draw_region; 

Rect empty_rect = (0, 0, 0, 0); 

RgnHandle empty_rgn; 

/« V 

void 

dcucontrolsO 

{ 

) 

/« .•! 

void 

dQ_menu(command) 
long command; 

{ 

int meniuid = HiWordfeommand); 
int item = LoWord(command); 
char item_name{32]; 
register int deLJiumber; 

5Witch(menu_id) 

( case APPLE-MENU: 

Ifdtem =« ABOUTJTEM) 

{ da_about_menuO; 

) 

else 

{ GetItem(GetMHandle<menu-id), item, item_name); 

OpenDeskAcc(item_name); 

da_number • ((WindowPeek)ProntWindow())->wlndowKind; 
if(da_number < 0 && -f^DESKOPEN 1) 
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) 

/» 


( EnableUeTn(QetMHanaie<EDlT_MENU), 0>; 

) 

) 

break; 

case FILE-MENU: 

daJUesdtem); 

break; 

case EDIT-MENU; 

SystemEditdtem - I); 
break; 

case PLASMA-MENU, 
swltchdtem) 

{ case INITIALIZE: 
break; 

case EXAMINE: 

MORE-PLOTS = TRUE; 
do-examineO; 
case INTERRUPT: 
break; 

case TRANSPARENT: 
toggleijeeO; 
break; 

case RESTART: 

RETURN = TRUE; 
break; 

) 

) 

HiliteMenu(O); 


void 

douabout-menuO 

( 

) 

/• 

void 

do-niesdtem) 


•V 


V 


Int item; 

{ register int de^jaumber « <(WindowPeek)ProntWlndow<»->windowKind; 
switchdtem) 

{ case PRINT: 

if(deL_number >= 0) 

( do-prlntO; 

lf(PrlntErr »= 0) 

{ handleL_error(PrlntErr); 

) 

} 

break; 
case NEW: 

it(da-number >= 0) 

{ do-ncwO; 

) 

break; 
case OPEN: 

U(da-number >= 0) 

{ doL.open(draw_windowt3], &draw-picture[3]); 
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) 

brsak; 
case CLOSE; 

if(da_number < 0) 

( CloseDeskActKda-number); 

if(--DESKOPEN <= 0> 

{ Disableltem(GetMHandle<E01T_Mf:NU), 0); 

) 

> 

break; 
case SAVE: 
case SAVEAS. 

if(da_number >= 0) 

{ dOL-save_asO; 

) 

break; 
case QUIT: 

finlshO; 

default: 

break; 

) 

) 

/. V 

void 

dOLjnouse^own(event) 

register EventRecord ‘event; 

{ register GrafPtr save_graf; 

Point b; 

WlndowPtr mouse. window; 
register int da.number; 

Int place-type » PlndWlndow(event->where, Amouse. window); 
register long getw » GetWRefCon(mouse-window); 

Boolean test; 
lf(lempty_rgn) 

{ cmpty_rgn = NewRgnO; 

RectRgn(empty_rgn, &empty_rect>; 

) 

swltch(place-type) 

{ case InDesk: 

break; 

case InMenuBar: 

dOLjnenu(MenuSelect(event->where»; 

break; 

case inSysWindow: 

SystemCllckCevent, mouse. window); 

SystemTaskO; 

da_number = ((WlndowPeek)FrontWlndowO)->wlndowKind; 
lf(d 2 L.number >° 0 AA --DESROPEN <= 0) 

( Dlsableltem(GetMHandle<EDITJVlENU), 0); 

} 

break; 

case inContent: 

lf{ mouse. window 1= FrontWlndowO) 
SelectWlndow(mouseL.wlndow); 

else 

doucontrols(mouse.wlndow, event->where); 
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) 

/• 


breaK; 

CAM inDrag: 

lf((test (getw < 0 && mouscL-window == ProntWindowO))) 

( INI • see_through()i 

save_rect » overlap-rect; 

DRAG = TRUE; 
lf( INI ) 

{ InvalRecU&saveu'ect); 

vaIid_scroll3(mouseL. window); 
vaIid_grow(mouse-wlndow); 

) 

) 

DragWindow<mouse-window, event->where, fcdragj'ect); 

If(test) 

{ INF = s«e-through(); 

if((IINI) && INF && (EmptyRgn( 

((WlndowPeeR)mouseL.window)->updateRgn))) 
{ InvalRectC&overlapLxect); 

UNDER_ONLY = TRUE; 

) 

if((IINF) M. (UNI) && EmptyRgn( 

((WindowPeek)mouse-wlndow)->updateRgn)) 
{ DRAG = FALSE; 

} 

) 

break; 
case InGrow: 

ifUest « (getw < 0 && mouse- window »» FrontWindowO)) 

{ INI s see_through(); 

GROW » TRUE; 
save-rect = overlap-rect; 

} 

grow_wlndow(mou3e_window, event->whcrc); 
lf(£mptyRgn(((WindowPeek)mouse_wlndow)->updateRgn)) 

{ GROW = FALSE; 

break; 

) 

INF = secL-throughO; 
break; 

case InGoAway: 

It (TrackGoAway(mousc_wlndow, event->where)) 

( HideWlndow(mouse-window); 

) 

break; 

default: 

break; 


V 


void 

grow_u pda te( window) 

register WindowPtr window; 
{ GrafPtr save_graf; 
GetPort(Asave_graf); 
SetPort(wlndow); 
Ifddraw-reglon) 



Page 98 


( draw_reglon = NewRgnO; 

) 

RectRgn(draw_region. &overlap_rect); 

SectRgn(draw_region, «WindowPeek)window>->updateRgn, draw_region); 
3etPort(save-^af>i 

) 

f V 

void 

dcuupdate(event) 

register EventRecord 'event; 

( GrafPtr save^af; 

register WindowPtr update_windorw > (WindowPtr)event->message; 
register long getw = GetWRefCon(update_window); 

GetPorUtrsavcL^raf); 

SetPort(update_window); 

BeginUpdate(ui)date.window); 

EraseRect(&update_window->portRect); 
DrawPicture{(PicHand]eXn)y_abs(getw)), fcbig_rect); 
UpdtControl(updateL.window, update_window->vlsRgn); 
DrawGrowlcon(updateL.window); 

EndUpdate<updateL.window); 

SetPort(save_graf); 


) 

/« *( 

void 

dOL.see_update<event) 

register EventRecord 'event; 

{ GrafPtr save_graf; 

register WindowPtr updateL_wlndow = (WindowPtr)event->message; 
register long getw ° GetWRefCon(update_window); 

Rect comiLJ'ect; 

GetPorK&saveL^raf); 

SetPort(updatCL.wlndow); 

BeginUpdate(update_window); 

EraseRect(&update^window->portRect); 
DrawPicture«PicHandleXmy-abs(getw)), kbig-rect); 
UpdtControKupdate. window, update_wlndow->visRgn); 
OrawGrowIcon(updateL_wlndow); 

EndUpdate(update-window); 

SetPort(save_graf); 
if(UNP s see_through(») 

{ lnvalRect(&overIap_rect); 
if(GROW) 

{ ValldRect(&save-rect); 

/'SectRectfJisave^ow, tioverlapurect, isaveLgrow);'/ 
InvaiRecti&savcLgrow); 

SectRecKtcsaveiJi, Goverlaii-rect, tcsaveLh); 
InvalRect(&savcLJi); 

SectRect(&sav«L.v, ftoverlaji-rect, &save_v); 
InvalRect(&saveL.v); 

) 

valld-scrons(update_ window); 
valld_grow<update- window); 

UNOER_ONLY - TRUE; 

if(£mptyRgn(((WlndowPeeK)updateL.window)->updateRgn)) 
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( UNOER-ONLY - 6ROW - FALSE; 

) 

) 

) 

/. V 

void 

under_update(event) 

raster EventRecord ‘event; 

( Point origin; 

GrafPtr sAvei^af; 

register WlndowPtr updatei-window - (WindowPtr)event->mess&ge: 
GetPDrt(&saveL 4 iraf); 

BeginUpdate(updateL.wlndow); 

/■EraseRcct(fcsaveL^ow);V 
SetPort(under-window); 
orlgin.h “ origin.v « 0; 

LocalToOlobaKAorlgln); 

SetPort(top-wlndow); 

GlobalToLocal(tiorigin); 

SetOrlgln(-origln.h, -origin.v); 

DrawPicture«PicHandle)(my-abs(GetWSefCon(under_wlndow))), Abig-rect); 
Set0rlgln(0, 0); 

DrawPlcture((PlcHandleXmy_abs(GetWRefCon(update-wlndow))),&big_rect); 

EndUpdate(updateL.vvrindow); 

SetPort(save_graf); 

UNDER-ONLY » DRAG > GROW = FALSE; 

) 

/. V 


void 

dauactivate<event) 

register EventRecord ‘event; 

( register WlndowPtr event-window » (WlndowPtrXevent->message>; 
GrafPtr saveLgraf; 
if(actlveFlag & event->modifiers) 

{ SetPort(event_window>; 

if(GetWRefCon(event-window) < 0) 

( INF B see_tbroughO; 

lnvalRect(&overlap_rect); 

valid-scroUs(event-window); 

) 

show-scrolls(event_window); 

) 

else 

( if(GetWRefCon(event-wlndow) < 0) 

( GetPort(AsaveL 4 raf); 

SetPort(event_wlndow); 

InvalRecU&overlaii-rect); 
valid-scrolls(event- window); 

SetPorUsaveLgraf); 

) 

hide-scroUs(event-window); 

) 

DrawGrowicon(event-window); 

/•valld-grow<event-wlndow);‘/ 

) 
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I* «/ 

void 

vftlld_rect(window, bod_rect) 

register WindowPtr window; 
register Rect *bad_rect; 

{ register ControlHandie control = ((WindowPeeK)window>->contrblList; 
register iong iabei; 

Rect vaLrect; 
while<control) 

( iabel •> GetCRefContcontrol); 

ifOabel == Uong)VjSCROLL I label Oong)H_SCROLL) 

{ SectRect(&(*control}->contrlRect, bad_rect, &val_rect); 
ValidRect(t>val_rect); 

) 

control = (•control)->nextControl; 

) 

/*defineLgrow<window, &val_rect); 

SectRect(&val_rect, bad_rect, &val_rect); 

ValidRect(&val_rect);V 

) 

/» . V 

void 

valid_scrolls(window) 

register WindowPtr window; 

{ register ControlHandie control = ((WindowPeeK)wlndow)->controlList; 
register long label; 
while(control) 

( label a GetCRefCon(control); 

Ifdabel (long)V_SCROU. I label (long)H_SCR01X) 

( ValidRect(gt<*control>->contrlRect); 

} 

control = (•control)->nextControl; 

) 

) 

I* V 

void 

show_scrolls(wlndow) 

register WindowPtr window; 

( register ControlHandie control ° ((WindowPeek)window)->controlList; 
register long label; 
while(control) 

( label > GetCRefCon(control); 

ifOabel 0ong)VJ9CROLL I label " Oong)H_SCROLL) 

{ HillteControKcontrol, 0); 

ValidRect(St(*control]^>contrlRect); 

) 

control “ (•control)->nextControl; 

) 

) 

/. V 

void 

hide_8crolls(window) 

register WindowPtr window; 

{ register ControlHandie control « ((WindowPeek)window)->controlList; 
register long label; 
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wbUe(control) 

( labal = GatCBafCon(control); 

ifUabel Uong)V^BOLL I label OonglHLSCSOLL) 

( HiliteControl(control, 255); 

Valld8cct(&(*control>->contrlRect); 

) 

control = (*control)->nextControl; 

} 

) 

/. •/ 

void 

moveLJ>ars(wlndow) 

register Windowptr window; 

{ 

register WlndowPeek peek = (WlndowPeek)window; 

register ControlHandle control = peek->controlLl3t; 

Int new.top = wlndow->portRect.top; 

Int newJeft ■ wlndow->portRect.left; 
int new_bottom = wlndow->portRect .bottom; 

Int newj'ight = wlndow->portRect. right; 
register long label; 

while (control) 

( label a GetCRefContcontrol); 
if (label (long)V-SCROLL) 

( HldcControKcontrol); 

save_v a (»control)->contrlRect; 

MoveControKcontrol, new_right - BAR-WIDTH, new_top -1); 
SlzeControKcontrol, 16, newJMttom - new_top - 15); 
ShowControKcontrol); 

ValldRect(fc((*control)->contrlRect)); 

) 

else if (label =a (long)M_SCROLL) 

{ HideControKcontroI); 

save_h = (*control)->contrlRect; 

MoveControKcontrol, new-left - 1, new_bottom - BAR-WIDTH); 
SlzeControKcontrol, new-rlght - new-left -13, 16); 
ShowControKcontrol); 

ValldRect(&((*control)->contrlRect)); 

) 

control • (*control)->nextControl; 

) 

) 

f* V 

void 

grow_wlndow(wlndow, mouse-point) 
register WlndowPtr window; 

Point mouse-point; 

{ long new-bounds; 

new-bounds ■ GrowWlndowfwlndow, mouse-point, Agrow-bounds); 
lf(new-J»ounds 1° 0) 

( invaLgrow(window); 

/*EraseRect(Asave-grow);*/ 

SlzeWlndow(window. LoWord(new-bounds), HiWord(new-bounds), TRUE); 
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move-bars(wlndow}; 

DrawGrowIcon<wlndow)i 

valld_grow<wlndow); 

) 

) 

/. V 

void 

deflne_grow<wlndow, grow_rect) 

WlndowPtr window; 
register Rect *grow_rect; 

( Rect port_rect; 

port-rect * wlndow->portRect; 

SetRect(grow_rect, port_rect .right - 16, port_rect. bottom - 16, 
port_rect. right, port-rect.bottom); 

) 

/» •/ 

void 

lnval_grow< window) 

WlndowPtr window; 

( Rect temp_rect; 

deflne_grow(wlndow, &temp-rect); 

InvalRect(&temp_rect); 
save_grow = temp_rect; 

} 

/» V 

void 

valid_grow< window) 

WlndowPtr window; 

( Rect temp-rect; 

define_grow<wlndow, &temp_rect); 

ValldRect(8( temp-rect); 

} 

/. */ 

void 

do_examlne<) 

( int response; 

static Int window-number - 0; 

Int selection; 

register WlndowPtr the_ window; 
if(my_ref == -999) 

( selection > window-number X 4; 

daL.open(draw-window(selectlon}, Sidraw-picturelselection] ); 

/•SelectWlndow<draw-window(selectlon]);V 

wlndow-number++; 

) 

whlle<MORE_PLOTS) 

( thcL-window » FrontWlndowO; 

hideL^rolls(theL.window); 
eraseL^row<the_wlndow); 
response = Alert(257, OL); 
show-scrolls(theL.window); 

OrawGrowIcon(the-wlndow); 

swltch(response) 

( case ENOUGH; 

MORE-PLOTS • FALSE; 
selection = window-number * 4; 
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dOLjopenWraw-windowfselectlon], Mraw_plcture(selecUon] ); 

/*S«lectWindow(draw-wlndcrw(Mlaction]);V 

da-window_care<); 

■window_nuinber++; 

break; 

case NOT-ENOUGH; 

selection « window_number X 4; 

/*SelectWindow<draw_wlndowlselection]);V 

do_read(draw_window(selecUon], 

&draw_picture[selecUonl, reply, my_ref); 
dcL.window_care<); 
window jiumber++; 
break; 

) 

) 

) 

/. 

void 

da-wlndow_care() 

( Boolean activate ° TRUE; 

EventRecord thia^vent; 
whlle<actlvate) 

{ EventAvalKeveryEvent, AthlSu^ent); 

IftthlSL^ent.what actlvateEvt) 

{ GetNextEventleveryEvcnt, atblSL^ent); 

da_actlvate<&thls-event); 

) 

else 

{ activate <= FALSE; 

} 

) 

} 

void 

erasc_grow( window) 

register WindowPtr window; 

{ Rect tempLj'ect, port-rect; 

port_rect = wlndow->portRect; 

SetRect(&tempj'ect, port_rect.right - 14, port_rect.bottom - 14, 
port-rect.right, port_rect.bottom); 

EraseRect(&temp_rect); 

) 

f ,, 

void 

toggle_see<) 

{ topL.window » ProntWlndowO; 

SetWRefConCtopuwlndow, -GetWRefCon(top_wlndow)}; 

IKUNl ■ INF s see^throu^O)) 

( lnvalRect<&overlap_rect); 
velid_scrolls(top- window); 
valld_grow<top_wlndow); 
lf(GetWRefCon(top-window) < 0) 

{ UNDER_ONLY - TRUE; 

) 

) 

) 
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/• V 

Boolean 
see_throughO 
{ GrafPtr save_graf; 

Boolean IN • TRUE^ 
top_wlndow = FrontWindowO; 

under_wlndow = (WlndowPtrX((WindowPeeKHop_wlndow)->nextWindo(w); 
GetPorU&save^af); 

( SetPort(under.wlndow); 

Loc 2 dToGlobal((Point *X&under_wlndow->portRect))i 
LocalToGlobal((Polnt “Xfc((under_wlndow->portRect).bottom))); 
SetPort(top-window); 

LocalToGlobal((Polnt *Xfctop_wlndow->portRect»i 
LocalToGlobaK(Point •Xfc((lopuwindow->porlRecl).botlom))); 
SectRect(&<top-window->portRect), &(under_wlndow->portRect}, 
tioverlaii-rect); 

lf(((long "iJcOverlapL-rectXO] == OL M ((long *)4overlap_rectXl] = OL) 
{ IN = FALSE; 

) 

GlobalToLocal((Point *X&top-window->portRect»; 

GlobalToLocal((Polnt *X&((top-window->portRectibottom))); 
GIobalToLocal((Polnt *X&overlap_rect)); 

GlobalToLocal((Poin t *X&over lap-r ect . bottom)); 
SetPort(under.window); 

GlobalToLocal((Point •XAunder_wlndow->portRect)); 
61obalToLocal((Point •X8i((undcr_wlndow->portRectibottom))); 
SetPort(saveL_graf); 

) 

return( IN ); 


) 



Page 105 


VI. advance.c 

This module contains the heart of the simulation program. The routines 
that carry out the “leap frog” advancement of the simulation particles in time and 
the routines that calculate and advance the field quantities are in this module. 
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*lncluae ‘piasma.h* 
extern double *x. “vx, “vy; 

extern double ael, epsi, rhoO, rho(NGlM], philNGlM], c{NGlM], eO, wO; 

extern int it, Ithl, irho, Irhos, iphi, ie, Ixvx, ivxvy, ifvx; 

extern double elapsed_Ume; 

extern double acratchlNGlM], fftUNGTWOl 

extern int ng, iw, 

extern double L, dx, dt; 

extern double e^NTHl], nms(NSPMl 

extern double *eseni[MMAX]; 

extern int mplotiMMAX]; 

extern double x_polnts[NGMAX^ 

extern WindowPtr draw_wlndow[4^ 

extern PicHandle draw-pictureCd]; 

extern int rho_file. rhoaJile, e_file, phi-file, vxvy_file, xvx_file, 
fv-file; 

void getJieldsO, transform_rhoO, rho-to-phiO, saveijaaodesO, doJnverseO; 
void electric-fieldO, reset-rhoO, accelerate, nearest-pointO; 
void p.xonserve(), e^xonserveO, advanceO; 

extern void field_initO, sincos20, inverscL^incos20, savei^aphOi 
extern void trace_plotO, showO; 

extern PILE *output_file; 

void 

get-fieldsUth) 

register int ith; 

{ static int ng2 = 0; 

sUUc double sm(NG2M], ks<|l[NG2Ml 
static max_min eLJimit^2]; 
static max_min rha-llmits(2l 
static max_mln phl_limits(2i 
double Li; 
register int i; 
if(lng2) 

( ng2 * ng/2; 

field-iniUsm, ksql, ng2); 

) 

rbo(0] rtao(ng]; 
rho(n^ « rho(0]; 

Iftirbo !■ 0 && Otta X irho =* 0)) 

{ save_grapta(rho-file, x-polnts, rho, rtaaJlmlts, ng, draw_window(0], 

g4raw_picture{oi, *\pCharge Density', TRUE, tretce_plot); 

) 

transform-rhoO; 
rbo-to-phiUth, ng2, ksqi, sm); 
save_modes(lth): 
dOLjnverseO; 

Iftirhos 1° 0 M Gth X irhos «■= 0)) 

( save_graph(rhos_file, x.polnts, rho, rhoulimits, ng, draw_window(0], 
&draw_plcture[0], *\pSmoothed Charge Density*, TRUE, trace-plot); 

) 

UGphi l> 0 U Gth X iphi == 0» 

( savcL^aphtphi-file, x_points, phi, phi-limits, ng, draw.windowU], 
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&draw_plcture{l], *\pPotentlal Phi", TRUE, trace-plot); 

) 

electrlc-fieldO; 

ifde l« 0 M. (1th X ie » 0)) 

( 3aveL4raph(e-ft)e, x_polnts, e, eJlmits, ng, draw_wlndow{2], 
&draw_plctur^2], “\pElectrlc Field*, TRUE, trace.plot); 

) 

resetJ'hoO; 

/• fprlntf(output_flle,*\nreset_rho*); 

show<x-polnts, rho, ng * 1);*/ 
ael c 1.; 

) 

/» .__v 

void 

transform-rhoO 
( register Int 1; 

double hdx; 

register double *rhoK ° rho; 

hdx <■ .5 * dx; 

for( 1 - 0; 1 < ng; !♦♦) 

{ rho(i] hdx; 
scratch[i] - 0.; 

) 

slncos2(rho, scratch, fftl, rhoK, scratch, ng); 
rhok[0] = 0.; 

) 

/« V 

void 

rhoutouphiUth, ng2, ksqi, sm) 
int ith, ng2; 
register double *ksql; 
double *sm; 

( double eses « 0.; 

register int i, index; 
register double *phik phi; 
register double *rhok ^ rho; 
phlk(0] = 0.; 
index > ng - 1; 

for ( 1 ■ 1; 1 < ng2; !♦♦, index--) 

( phik[l] « ksqiCl - I] * rhok[i]; 

phikdndex] » ksqill - 1] * rhoklindex]; 
eses rhokUl * pbikll] * rhokllndex] * phikdndex]; 
rhok(l] smll - IJ; 
rhokdndex] *= smd - 1); 

) 

phik[ng2] ° ksqilng2 - 1] * rhok{ng2]; 

eseCith] ° (2. * eses ♦ rhok[ng2] * phiklng2]) / (2. * L); 

rhok{ng2] sm(ng2 - IJ; 

) 

/« V 

void 

save-modesdth) 

register int ith; 

( register int 1 » 0; 

register int index - ng; 
double temp; 
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regisUr double *phik ■ phi; 
register double *rhok = rho; 
whlle<(inplotti] l« 0) M a »= MMAX)) 

{ temp ■ (rhoklO * phlkti] + rhokUndex] * phikllndexP / L; 
ifG °o index) 

{ temp •= .25; 

) 

esemO^^Iith] = temp; 
index--; 

) 

/• •/ 

void 

dQ_inverse<) 

( register double *phik ° phi; 
register double *rhok = rho; 
double Li « 1. / L; 
register int i; 
for( i = 0; i < ng; !♦♦) 

{ rhokli] “= Li; 
phikli] •= U; 

) 

inverse^incos2(pbik, rhok, fftl, phi, rho, ng); 
phiing] > phi(Ol; 
rho(ng] = rhoCOl; 

/• fprinti(output-file, *\ndouinverse*); 

show(x_points, rho, ng ♦ l);V 

) 

/. V 

void 

electrlcJieldO 
{ double eOt, hdxl, dxi; 

register int i; 

eOt B eo * cos(wO * elapsed-tlme); 
swltchOw) 

( case ZERO-ORDER: 
case MOMENTUM: 
hdxi = .5 / dx; 
lor (1 - 1; 1 < ng; !♦♦) 

( e{ij » (phl[i-l] - phi[i + 1]) • hdxl ♦ eOt; 

) 

elO] « (phKng - 13 - phi[l]) • hdxi ♦ eOt; 
e{ngj = e{Oj; 
break; 
case ENERGY: 

dxi B I. / dx; 

for (1 » 0; 1 < ng; i++) 

{ elll = (phlll] - phltl *11) • dxi ♦ eOt; 

} 

e(ng] = e{0]; 
break; 

) 

) 

/. V 

void 

reset_rho() 
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( register Int 1; 

for ( 1 B Oj i < ng; !♦♦) 

{ rhoii] B rhoO; 

) 

rho[ng] ■ 0.; 

) 

/« „.„V 

void 

acceleratedl, iu, q, m, t, p, ke> 
double q, m, 
register double *p, *ke; 
register Int 11, iui 
{ register double ‘ace; 
double dxdt = dx / dt; 
double ae B (q / m) * (dt / dxdt); 
double temp; 
register int 1; 
lf(NEED-R0TAT10N) 

( ae •» .5; 

) 

if(ae Is ael) 

( temp s ae / ael; 

for( 1 = 0; 1 <= ng; !♦♦) 

{ atU •= temp; 

) 

ael s ae; 

} 

swltchOw) 

( case 2 ESOlOROER; 

nearest_polnt(il, iu, dxdt, m, ke, p>; 
break; 

case MOMENTUM; 

puconserveOl, iu, dxdt, m, t, ke, p); 
break; 
case ENERGY: 

eL.conserve(il, iu, dxdt, m, ke, p); 
break; 

) 

} 

/» V 

void 

nearest^pointOl, iu, dxdt, m, ke, p) 
register int U, iu; 
double dxdt, m; 
register double *ke, *p; 

{ double vis ■ 0.; 
double v2s • 0.; 
register double *a ° e; 
double vO, vn; 
intl.j; 

for( 1 • 11; 1 < iu; i+*) 

{ J - (xtll ♦ .5); 
vO » vxli]; 
vn = vO ♦ a£jt 
vis ♦» vn; 
v2s ♦= vn • vO; 
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vxUl » vn; 

) 

•p ♦= m • vis ■ dxdt; 

•ke ♦- .5 • m • v2s • dxdt “ dxdt; 

) 

/• V 

void 

P-cons«rve<il, lu, dxdt, m, t, ke, p) 
register Int 11, lu; 
double dxdt, m, ti 
register double *p, *ke; 

( register int 1, J; 

double vO, vn, vis, v2s, s, aa, vxx, vyy; 
register double *a « e; 

If(NEEBLROTATION) 

{ s - 2. * t / ( 1. ♦ t • t); 
v2s » 0.; 

fori 1 » il; 1 < lu; !♦♦) 

{ J = xIH 

aa = aljj ♦ (x(l] - J) ■ (a[J ♦ 1] - alj]); 
vyy = vydj; 

vxx = vx{i] - t • vyy ♦ aa; 

vyy ♦= s • vxx; 

vxx -= t “ vyy; 

v2s ♦= vxx • vxx ♦ vyy • vyy; 

vxtl] = vxx ♦ aa; 

vy(i) = vyy; 

) 

•ke ♦» .5 “ m • v2s • dxdt • dxdt; 

) 

else 

( vis » 0.; 
v2s = 0.; 

for( i s 11; 1 < iu; !♦♦) 

( J " xUl 
vO >» vx[lj; 

vn ■ vO ♦ a[J] ♦ CxCl] - J) * (aU ♦ 1] - alj]); 
vis vn; 
v2s vO • vn; 
vx(l] = vn; 

) 

•p m ■ vis * dxdt; 

•ke .5 • m • v2s • dxdt • dxdt; 

) 

/» •/ 

void 

e_conserve(il, iu, dxdt, m, ke, p) 
register Int 11, lu; 
double dxdt, m; 
register double *p, *ke; 

( register int 1, J; 

double vO, vn, vis, v2s; 
register double *a ■ e; 
vis = v2s = 0.; 
fori 1 ■ 11; 1 < iu; !♦♦) 



/ 
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t J “ x[H 

vO » vx(i]i 

vn » vO ♦ aljl; 
vis ♦» vn; 
v2s ♦» vn • vO; 
vxlll » vn; 

) 

•p ♦=» m * vis ■ dxdt; 

•ke ♦« .5 • m • v2s • dxdt • dxdt; 

) 

/. 

void 

advancedl, iu, q) 

register int II, iu; 
double q; 

{ register int 1, J; 

double qdx « q / dx; 
double xn » ng; 
double drho; 
switchUw) 

( case ZESO_OKOER: 

tori i - 11; i < iu; !+♦) 

{ xlil ♦= vxtl]; 
if(x{i) < 0.) 

{ x[l] ♦= xn; 

) 

if(x(l) >= xn) 

( x[l} -• xn; 

) 

J « xtll ♦ .5; 
rhotj] ♦“ qdx; 

) 

break; 

case MOMENTUM: 
case ENERGY: 

for( i = il; 1 < iu; !♦+) 

{ x[i] vxlll; 

U(xll] < 0.) 

{ xli] ♦" xn; 

) 

if(x(i] >= xn) 

{ xli] xn; 

) 

i - klO; 

drho = qdx ■ (xW - J); 
rhoUl ♦“ qdx - drho; 
rho(j ♦ 1) ♦= drho; 

) 

break; 

) 

) 


V 


X 
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VII, fourier.c 

This module consists of the fast Fourier transform routines utilized in 
advance.c. 



^include *|>lasma.b* 
double sine«{lS], slnepdS]; 

void fULjsine^tableO. fourlO. sinco520. inverscL^incos2(); 
void 

fill_sine_tableO 
{ double t a atand.); 
register int i^ 
sines(0] a 1.; 
for(i a 1; 1 < 15; !♦♦) 

{ Sinead] a sin(t); 

sinep(i] a -2. • powtsinesli], 2.); 

t •« .5; 

) 

) 

void 

fourKdata, nn, isign) 

register double *data; 
int nn, isign; 

( register int i; 

register int j ■ 0; 

int n a 2 • nn; 

register int mmax, m, istep; 

int sincnt; 

double wtemp, wr, wpr, wpi, wi, tempr, tempi; 
ford a 0; i < n; i ♦= 2) 

( if(J > i) 

( tempr a dataUli 

tempi a dataU * 1]; 
dautj] a datali]; 
dataCi 1] a daiali * 1]; 
data(i] a tempr; 
data[i + 1] a tempi; 

) 

m a nn; 

while< (m >a 2) M ((J ^ 1) > m)) 

{ J -= m; 

m /a 2; 

) 

J «n; 

) 

mmax a 2; 
wpi a 0.; 
sincnt a 0; 
whilednmax < n) 

( istep a 2 * mmax; 
wpr a sinep(sincnt% 
wr a 1,J 
wi a 0.; 

fordn a 0; m < mmax; m '»a2) 

{ ford a m; 1 < n; i ♦= istep) 

( J a i « mmax; 

tempr = wr • data(J] - wi • dataU ♦ 1]; 
tempi a wr • datalj ♦ 1] ♦ wl • data[J]t 
dat^] a data[i] - tempr; 
data(J 1] a datali 1] - tempi; 
datali] tempr; 



datt^i ♦ 1 ] tempi: 


) 

wtemp ■> wr; 

wr ♦* wr “ wpr - wi • wpi; 
wl ♦= wl • wpr ♦ wtemp • wpi; 

) 

mmax •= 2; 

wpl « (Isign > 0? stneslslncnt] ; -slneslslncnt]): 
slncnt««: 

) 

) 

void 

slncos2(datal, data2, fft, datakl, datak2, n) 
double *datal, *data2; 
register double *datakl, *datak2; 
register double (*fftX2): 
register Int n; 

{ double rp, ip, rm, im; 
register int J, index; 
foHj « 0; J < n; J++) 

( fftIJIO] datalCj]; 

ffUjll] » dato2U^ 

) 

(ourlCfft, n, 1); 

datakUO] ° 2. * ffUOlO]; 

datak2[0] 2. * flUOll]; 

for(J = 1, index “ n - 1; J < index; J++, index—) 

( rp - tfUjlOj; 

ip = tttum 

rm « ffUindexlO]; 

Im - tfUindexIl]; 
datakUj] * (rp ♦ rm^, 
datak2[J] = Up ■» im); 
datakl[index] » Up - im); 
datak2[index] « (rm - rp); 

) 

dataklU] ° 2. * f/t(Jl0]; 
datak2[J] « 2. * ffUjIl]; 

) 

void 

inverse_slncos2(datakl, datak2. itt, datal, data2, n) 
double *datal, *data2i 
register double *datakl, *datak2; 
register double (*fftX2l 
register int n; 

( double ca, sb, cb, sa; 
register Int J, index; 

(ftlOlO] - datakUO]; 
ffUOlI] > daUk2[0}; 

lor(4 = 1, Index = n - 1; j < index; >♦, index—) 

( ca o datakUj]; 
cb ■ datak2[J]; 
sa > datakUindex]; 
sb > datak2[lndexi 
(fUjIO] - ca - sb; 
fftUIll o cb * sa; 



f/t{lndexl 0 } a ca « 
^ f/Uind»xIl} = cb - 

f/ttjlo] - datakl[j]; 
MtlJIll = dat&K2{j]; 
fourJ(//t, n, -l>, 
forij = 0; J < n; >♦) 

{ datallj] = fftUIOi 

daU2[j] a «tOIH 
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VIII. plajsma.h 

This is the header file that contains definitions of various quantities used in 
the other modules. Many of these definitions are in the form of include 
statements for still other header files, mainly 


Macintosh system and i/o definitions. 
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•inciuae "matti.h" 

^include ‘stdio.h* 

*include “unlx.h' 

^include ‘sanc.h* 

*includs *QuickDraw.h* 

•include ‘MemoryMgr.h’ 

•include ‘MacT/^.h* 

•include “WlndowMgr.h* 

•include “ControlMgr.h* 

•include “EventMgr.h* 

•include ‘DesKMgr.h* 

•include 'MenuMgr.h* 

•include *ToolboxUtll.h* 

•include ‘DlalogMgr.h* 

•include ‘ResourceMgr.h* 

•include ‘FontMgr.h* 

•include ■TextEdlt.h" 

•include “PrlntMgr.h" 

•include •PlleMgr.h* 

•Include ‘StdPilePkg.h* 

•define my_abs(x) «x) > 07 (x> ; -(x» 
•define my_max<x,y) ((x) > (y)7 (x) : (y» 

•define NTH 500 
•define MMAX 16 
•define NSPM 3 
•define NGMAX 2S6 
•define NGIM (NGMAX * I) 

•define NG2M (NGMAX / 2) 

•define NGTWO (2 • NGMAX) 

•define NPAR 512 
•define NTHl (NTH ♦ 1) 

•define NTH2 (NTH ♦ 2) 

•define NSPMl (NSPM « I) 

•define ZERO.ORDER 1 
•define MOMENTUM 2 
•define ENERGY 3 

•define ORDERED (vt2 l>= 0.) 

•define MAGNETIZED (wc l> 0.) 

•define ANOTHER (nig !■ 0) 

•define RANDOM (vtl 1° 0.) 

•define PIRST_GROUP (11 »» 0) 

•define NEED-ROTATION (t 0.) 

•define V-SCROLL 256 
•define H_9CROLL 257 
•define BAR.WIDTH 15 
•define GRAPH-WINDOW 256 

•define APPLE.MENU 1 
•define PILE-MENU 256 
•define PLASMA-MENU 257 
•define EDIT-MENU 258 
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•dofin* ABOUT-ITEM 1 

•define NEW 1 
•define OPEN 2 
•define CLOSE 3 
•define SAVE 4 
•define SAVEAS 5 
•define REVERT 6 
•define PRINT 8 
•define QUIT 10 

•define ENDALERT 256 

•define INITIALIZE 1 
•define EXAMINE 2 
•define INTERRUPT 3 
•define TRANSPARENT 4 
•define RESTART 5 

•define DRVR Ox44525652L /» The string *DRVR* as a long */ 

•define mk_long(x) ("{Oong •) A (x))) 

typedef struct(double fmax, 

fmin;} hiax-mln; 
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9.0 PROBLEMS ENCOUNTERED DURING REPORTING PERIOD, R.D. 

ESTES, PI 

None 

10.0 ACTIVITY PLANNED FOR NEXT REPORTING PERIOD, R.D. ESTES, 

PI 

The work of this task is proceeding on schedule. In the coming period the 
software development will continue with the modification of the plasma simulation 
program to more closely approximate the physics of our hollow cathode plasmas. 
This means that collisional effects will have to be added and that the boundary 
conditions will have to be modified. Initial computer experiments will be made 
and their results compared with actual laboratory results. Parameters will be 
varied to determine which ones are most important. In addition to continuing the 
line of development already begun, we will be investigating the use of another 
computer code for applicability to our problems. MASK, a two-dimensional, fully 
electromagnetic, relativistic plasma simulation code developed at MIT and 
improved by Scientific Applications International Corp. (SAIC) is available for our 
use. The relativistic features are clearly not needed and would carry a substantial 
overhead, but it may be possible to modify the code easily to run in a 
nonrelativistic mode. A CONVEX computer is now available to us in the Radio 
and Geoastronomy Division of the Harvard-Smithsonian Center for Astrophysics. 
The vectorizing capabilities of this machine may make simulations with larger 
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